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ABSTRACT 


A digital computer program, originally developed by George E. Fink, 
capable of determining the in-and/or out-of-plane vibration frequencies 
of a single plane piping system, using the basic transfer method, is 
modified and augmented by the writer. The program is modified to use the 
writer's method to reduce the amount of calculation and is augmented to 
use Marguerre and Uhrig's modifying frequency determinant method and 
Pestel and Mahrenholtz's remainder method for higher natural frequencies. 

An analysis is made of difficulties encountered with the original 
transfer method and utilizing the modified methods. 

The program accepts branched systems and non-stiff intermediate 
supports. A discussion and an explanation of how to use the program are 


included. 
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CHAPTER I 


INTRODUCTION 


1-1 General remarks. 

Today, in increasing measure, natural frequencies of a system capa- 
ble of linear vibration are being determined with the help of the method 
of transfer matrices. This method is well explained in the recent book 
by E. Pestel and Leckie [3].* 

This method which has been known for more than ten years, proceeds 
with particular simplicity especially with the use of digital computer. 
However, it meets with important numerical difficulties. 

Numerical difficulties were observed by George E. Fink in his U. S. 
Naval Postgraduate School Master's thesis "Vibration analysis of piping 
ae He used the straight-forward method of transfer matrices 
for calculation of the natural frequencies of vibrating piping system$. 

The purpose of this paper is to deal with numerical difficulties 
and with the remedies which have been developed and also to investigate 


reduction in calculation. 


1-2 Scope of work. 

Numerical difficulties are encountered at higher natural frequencies 
and for systems having very stiff exterior springs. In examining these 
difficulties the paper by Marguerre and Uhrig was encountered [2]. This 
tends to explain difficulties both with the original transfer method and 
with the Pestel and Mahrenholtz's modified methods [1]. Although, 


Marguerre and Uhrig suggest several different methods for dealing with 


*Numbers in brackets refer to the bibliography, page 37. 


the difficulties, those which appear to be certain to eliminate the 
difficulties are so far outside the spirit of the transfer method that 
the advantage of this method would be lost. Therefore, it was decided 
to proceed with the transfer method, to get as much good from it as one 
could. 

The writer set out to investigate the remedies for failure at high 
frequency. The digital computer program "VIPIPE" presented in the Ap- 
pendices was used to determine the natural frequencies of in-and out-of- 
plane vibration of a planar piping system based on the transfer method. 
This program was originally developed by Fink [4], but has been modified 


and augmented by the writer. 


1-3 Assumptions and limitations. 

1. System is conservative. 

2. Distributed mass or/and lumped mass model is used. For the lump- 
ed mass model used in curved section (or in straight section) sufficient 
subsectioning is carried out to give acceptable answer. 

3. The transfer method surely fails if the system has rigid exterior 
elements. The hangers are introduced as a support by a point matrix with 
appreciable linear and/or torsional compliance. 

4. Program VIPIPE is limited to a system of a maximum of fifty com- 
ponents, twenty branches, and twenty hangers. It can only determine in- 
and/or out-of-plane natural frequencies of planar system. 

>. Boundary conditions should be such that half of the state quan- 


tities are vanishing.* 


*Refer to Appendix B. 


1-4 Notation. 


[ ] matrix. 
j column matrix. 

> HA cartesian right handed coordinate system. 
u,V,W displacement in the XYZ direction respectively. 
svv rotation about XYZ axis respectively. 

M bending moment. 

N normal force. 

T torque. 

V shear force. 

Z. state vector at location i. 

[Z,] state matrix at location i. 

[U,,J transfer matrix from i to i+ 1. 

A(w) frequency determinant. 

Am(W) modified determinant .* 
(N] non-vanishing state matrix of Zo 

Uj element of [U]. 

(s] spring matrix. 

[S] purged frequency determinant in the P-M method.* 
[I] unit matrix. 

[oO] null matrix. 

d displacement vector. 

{d] displacement matrix. 

P internal force vector. 

[p] internal force matrix. 


*In the text, the following abbreviations are used. M-D method stands 
for Modified Delta method and P-M method for Pestel and Mahrenholtz's 


remainder method. These will be explained in what follows. 
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[c] coordinate transfer matrix. 


[RI] y (R } J (R3} square sub-matrices. 

Vo column vector of every elements value of l. (cf. r 16 ) 
R(wW) remainder. | 

X correction factor vector. 


P» Py » Pp assumed state quantities at the starting boundary. 


X» X, » Xo correction factors in the P-M method. 


D, » Do purging factors in the M-D method. 
AX eigenvalue. 
x eigenvector. 








CHAPTER II 


TRANSFER METHODS AND NUMERICAL DIFFICULTIES 


2-l The state vector and transfer matrix. 
A. State vector. 
The state of vibration at location i, specified by state quan- 


tities, is described by the state vector Z This vector is a column 


i 
matrix whose elements completely describe the instaneous displacements, 
both rectilinear and angular, of that point from its quiescent position, 
as well as the corresponding rectilinear and angular forces in the member 
at the same point and at the same instant, i.e., those forces, if a cut- 
ting plane were passed through the point at the instant of interest, would 
have to be applied to each cut face to prevent relative motion between 
them. 
In the case of planar piping system, the state vector has six ele- 
ments, three pertaining to displacement, and three pertaining to force. 
The state vectors, for the in- and out-of-plane cases respectively, 
as used in VIPIPE, are: 
u displacement in X direction. 
Vv displacement in Y direction. 
1 slope in X-Y plane with respect to X asix. 
=|M; ‘moment about Z axis. (2-1-1) 
Vy shear in Y STP e Sei oem 


|N normal force in X direction. 


§ twist about X axis. 
T torque about X axis. 
-w (negative of) deflection in Z direction. 
Z.= slope in X-Y plane about X axis. ee) 
My moment about Y axis. 
Vy shear in Z direction. 
Thus, the state vector in this application is of order 6xl, while the 
system state matrix* is 6x3 and the transfer matrix of a section is 6x6; 
however, in developing theory we will use the more general case 2rxl, 2rxr, 


and 2rx2r respectively. A local XYZ right handed cartesian coordinate 


system and the elements of the state vector are arranged as in Fig. 2-1-1. 





Central axis 


Fig. 2-1-1 


B. Transfer matrix and transfer method. 


From the state vector at location i (node i), the state 


*The concept of "state matrix'' is introduced in Section 2-3. 
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vector at location i +1 is given by 


Zee, eae 2, Cl 0,1,2,...,0) (2-1-4) 


Here, matrix (Visi) relating the state vectors of node i and node i + 1 
is termed a transfer matrix. 

We number the node i + 1 immediately following the node i from the 
left end, the beginning of the main system. From the recurrence relation 


2-1-4, one easily obtains 


= (u) . Z, (2-1-5) 
which is a linear relationship between the state quantities at location 
O and n of the system (boundaries of the main system). We call this 
method the transfer method. 

Since of the 2r state quantities at each boundary r quantities are 
vanishing, only r are established by the boundary condition at location 
0, and r homogeneous equations are established by the boundary condition 
at location n. The natural frequencies of the system are given by the 
zeros of the determinant of the corresponding r homogeneous equations. 
Since the application of interest is to the determination of the natural 
frequencies of piping systems, it would be instructive to illustrate the 
origin and nature of frequency dependence in transfer matrices. For thjs 
end, let us construct a point matrix for the in-plane case connecting 7a 
with z* Ge. and Z; are state vectors of right side and left side, 


respectively, of a point mass m, at node i). 


i 
To begin we assume that all of the elements of the state vectors 


vary with time, each having the factor cos WL . 
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Noting that the deflection, slope, and moment* are continuous across 


the concentrated mass mi» so that 


L. 
= M: 


KR 
A L 


‘e 
uk =u- ; ve = vo ’ ute Ui ’ M ; 


From the free body diagram of Fig. 2-1-2A, B, the following relations 


are established. 


u(t) = u.cosWt 7 = -ul) cosWt (2-1-6) 
v(t) = v.cosy)t ot = -v\) coswWt (2-1-7) 
ve ° cos W t-V, cos lL) t = - TT: dv tt (2-1-8) 
NP + cos W) t-N; -cos|) t = - mM: us) (2-1-9) 


and substituting Eqn. 2-1-6 into Eqn. 2-1-8 and Eqn. 2-1-7 into Eqn. 


2-1-9 we have 


y. Vi +ovw (2-1-10) 
KR oo yk a 
Ne = N; + muy (2-1-11) 


Now, the relation of the state vectors of the right and left side of 


the point mass m in matrix form is: 


zi = (Up) - 2) 
iee. 5 

ees Sea wT! T u 

Vv 0 1 0 0 0 0 Vv 

v oo Ter oa ve 
M, 0 0 of Tr a my M3 
Vy 6“. oF “1% Vy 
Hie ay OOCOCOOtC CD N Je 


[up] is a point matrix for in-plane vibration. It neglects any rotational 


inertia that the point mass may have; however, rotational inertia could 


*The moment is continuous only if there are no point elements having finite 
rotational inertia. 
12 





be accounted for by including anothanvinonzero: term iin the ihatrix. 





L eo 3 
R mw U 
IM [vi pa 
a — 
|v; [vs My Ni Nj Ni 
mV 
Fig. 2-1-2A Fig. 2-1-2B 


The derivation of transfer matrices is treated in great detail in refer- 
ence (3) and transfer matrices for various components of planar piping 


system are incorporated in VIPIPE. 


2-2 Failure of transfer method. 

For higher natural frequencies and for systems having very stiff 
exterior springs numerical difficulties are encountered. In these cases, 
the numerical value of the frequency determinant is given as the differ- 
ence between large numbers. So practically the natural frequencies can 
no longer be determined. 

From a mathematical point of view, in case the transfer matrix [U] 
has a single dominant real eigenvalue, the columns of [u}* become more 
and more parallel as k increases, and approach the first eigenvector x) 


[6] [2] . In this case the transfer matrix degenerates. 


k “ ; 
[U] My = he XG Sa 


This could be the case for the lumped mass system as the number of lumped 
masses increases. This was observed in the straight section lumped mass 


“system, but for the curved lumped mass system this was not observed. 
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As the frequency increases (elements of some transfer matrices are 
functions of frequency), the same phenomenon was observed in the distri- 


buted mass straight section transfer matrix. For example, for the sample 


problem of Appendix E-2-1, 
At frequency of 20 rad/sec, 
the largest eigenvalue: A, = .11484022x107 


second largest eigenvalue: Darks = .99980414 + j0.01979160 


transfer matrix [U} is: 





99980413 —", 90535567\> 


— —— “t —_ as 22549093  . = 
ey 47565399" | .30070198.° | .25105587\'_ | .25986435\° . 2 
562196182 | .47565399\" | .30840073! | . 25105589. a 


x? means 1037, etc. 














In this case, that is, at this frequency, we can see that there is no 
phenomenon of the columns becoming parallel. As the frequency increases 
this single dominant real eigenvalue keeps increasing, and at frequency 
of 1700 rad/sec. 

the largest eigenvalue: A,= .5937220 x ig’? 

second largest eigenvalue: De d3 = -.11136847 + j.99377918 


transfer matrix [U] is: 
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-. 11136847 0 


53482189 > 


i 0 
| 148430788 | .13191200\14 : 
(0 ——_| -tes2ei4s\'5) re284e700 ,14843074\'°| .13191200\""| 0 J 
——* 16701804 » | 
0 


11136847 





0 





Normalized values of each column of [U) are: 


oni 
.9114641\"7 | .91146407 "7 | .9114641\77 0 
(aaa ea gE 
.1125225\7+ | .1125225\! | o 


ont 


0 / 0 i 






-8100280\"° 









We can see several of the columns are almost parallel to each other. From 
this, it is easy to see the numerical difficulties we should encounter. 
This depends on the nature of the transfer matrix. For systems made of 
various sections, multiplications of different transfer matrices may result 


in the difficulties taking care of themselves. 


2-3 Reduction in calculation in matrix multiplication. 


In equation 2-1-5, [U} = [U,]- [U,] ---(U,]. [Uj] » so we 
should perform matrix multiplication of 2rx2r matrix times 2rx2r matrix 


successively. 


Now using the fact of the boundary condition, introduce Z, as 


2rxr matrix called the "state matrix.'' This is formed from the state 
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vector as follows, all elements of the state matrix are zero except that 
in its ith column, unity appears in the same position as the ith general- 
ly nonzero term reading down the state vector. Then we can multiply 
2rx2r matrices by 2rxr matrices and thus can reduce the calculation 
considerably. 


For example, in the in-plane case, fora built-in left end 


a 0 0 0 
0 0 0 
0 0 0 
1 0 0 
0 1 0 
0 0 1 


The frequency determinant is obtained in the usual way from the boundary 
conditions at the right end. 

We use the term "Delta method" to refer to this method which was 
discovered by the writer during the course of this study. This pro- 
cedure is used in what follows instead of the usual procedure of the 


transfer method. 


2-4 Pastel and Mahrenholtz's remainder method. 
Instead of the frequency determinant, Pestel and Mahrenholtz intro- 


duced a remainder [1],[3] . This method, presenting the system state 





matrix as 2rxr matrix, reduces the amount of calculation considerably 
compared to the original transfer method, and claims greater accuracy at 
higher frequencies because the remainder, which is:the difference between 
small numbers, is itself small. This will be explained later. | 
A. Remainder. 
For the non-vanishing state quantities of the state vector Zo 
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a? 


which depends on frequency \) we use estimated values which approximate 
the true values, following which the Z. (i = 1, 2, 3, ...,n) are deter- 
mined. 

In the state vector Z,, since for eigenvibration only the ratios of 
state quantities are of interest, one of the nonvanishing State quanti- 
ties is assigned the value of l. 


For example, for the in-plane case, at certain frequency, a canti- 


levered pipe which is free at the left end and fixed at the right end has; 


+ X, (2-4-1') 


0 
(in what follows, superscripts 0 and 1 denote the round of calculation: 
i.e., super-script 0 means first round of calculation with approximated 
state quantities in Z, , and superscript 1 means second round of calcula- 
tion with improved state quantities in Z, ). In equation 2-4-1', we 
assigned the value of 1 to the non-vanishing state quantity in the first 
row of Z,, and the quantities of PY ; P, are the previously chosen 
approximations and os ; X° are correction factors which are unknowns. 
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The state vector Tip of (2-4-1') in other matrix form is: 


Z “ier Om oO} . [1 
Pl a 
rit Xs 
0 oOo 0 (2=4-1'") 
0 oOo 0 
Oo © 20 


ib 


For the actual matrix multiplication Zr, of (2-4-1'') is used. After a 

round of calculation, the approximated state quantities PY P. will 

be improved in the form 
Pee & Preche 


> “hh + (2-4-2) 


I 
} 
| ° 
P, 2 
where Ke and Kc are corrections obtained by this process; the 

difference between xe and x? ; xe and Ne will be shown later. 


Then, the improved state vector zy, will be: 
P, (2-4-3) 


0 
When we start first, usually P, 2 PS are set equal to l, in the 
absence of better information. From the same method as (2-1-5) using 
Z, of (2-4-1''). 
a ee ee oe eS 
From the boundary condition at location n, Le@., from the three vanishing 
elements of Z.. , we get te 
0 = A(W)-(N]* X= [S]) Xx (224-4 9 
[(N]} : nonvanishing sub-matrix of [Z,]. 
X : correction factor vector. 
A()) : frequency determinant. 
[s} : A@)-{n) 
For the cantilevered pipe of the example in the in-plane case, the 
frequency determinant A(W) , see Ref. [3] , consists of the top left 
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3x3 submatrix of the transfer matrix [UJ . 
Equation 2-4-4 may be written: 


0|= Uy Ui2 


Ui3 | e Q @ |e | 
4 Uri = Uge Ub P | 69 ¢ 
O} [Ur Ue Us} | Poo 1] L% 


Uy T Uia Pp 7 Us Pp. Uj Uys : 
Uy + U.P? + U23P. Use OU x, 


U3, t UP” t Ugs PL Uz. 33 X2 


Sai S22 $23 ee 
L $3 S532 533 Me 

In equation 2-4-5, the first column of (s] no longer contains large 
numbers, but the second and third columns remain unchanged. We call 


this "purging" the first column. From equation 2-4-5 we may write, 


Ver l= 1s. Sia | G 
5 Se aS G (2-4-6) 
23; S30 S33 


Equation 2-4-6 has first, three equations with two unknowns (two cor- 
rection factors), and second it has the purged column to the left and 
unpurged columns as coefficient matrix. 


Equation 2-4-6 is the same as; 


~~ Uy a Ui U3 fe P, 7h X, 
Uy) Uz. Up P, By (2-4-7) 
U3\ Usa U3 


From 2-4-7 we choose any two equations. If we choose the first two, 


then by Cramer's rule; 


We) 














—|Un Ups =e clalg; 
Uz, U3 0 A 0 é U22 Vat eB 
SS ae oe = eee 
Uin2 Urs | Uj2 ~ 
Ura asl U2. Up 


The A 's are defined by the preceeding expressions. If W=W, , 
where Wp is a natural frequency of the system, then the third equa- 
tion will be satisfied if we substitute these values for Xx, and K, ; 
Otherwise, substituting the above value for a , we get, from the 
third equation, a somewhat different value for XS ; we denote this by 


the symbol Y > 


Wee (2-4-9) 





We define the remainder R(0) to be the difference between these 


values. 


Ys) ° fo) a U3, A U3 
R(o)=X,-Y, = ra m. — — (224510) 





As we see from equation 2-4-10, for fixed \ , the remainder is indepen- 
dent of the assumed state quantities P and Pp, - If the remainder is 
zero, W = Wp since all the prescribed boundary conditions will be 
satisfied simultaneously. 

For the next round of caléulation for the same frequency, form the 


arithmatic average of the two values X and Y for the correction factor 
RR yore ss 

x= x, X= (+ 3) 
Then, new starting values are: 
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41 
| & o (as Us2 Ao 
ee ee — ae ai uk Wer 
Pa Pa 2 a a ea, ke ¢ 
which are also independent of Py and PS : A repeated round of cal- 
culation gives 
ray) 
no Ke gee we 
e-— = ot he 1 
ie ee 


el 
therefore, 
ee RG) 

Now, we see that the iteration is completed after the first round of 
calculation and the remainder is no longer the difference between very 
large numbers, because the remainder is of the same order of magnitude 
as the corrections xe, and Y, - During the iteration seeking the natural 
frequencies, we calculate R-(W)) and P; ( Wi) 's (i=1,2), and increase 
frequency by an amount At) to get Ido = Wt AW . Then get R (Wz) using 
Pe (Wz) =P. (1) (i = 1,2). The above example is for the case r = 3. In 
other cases for r >1, the spirit is the same (in case r = 1, we cannot 
use the P-M method). 

B. Various ways to get remainder. 

There are many ways to calculate a remainder. Each one has a 
different value at the same frequency (WEP), but passes through zero 
at W=Wp > 

Some fail and some give good results (see section 2-4-D). 
The possible ways arise from the following choices; 


21 


ie Which one of the nonvanishing state quantities will be 


given the value of 1 in the state vector Zo: 

2. Among the r homogeneous equations, as in equation 2-4-5, 
which r-l equations will be chosen to get x" : x; oe! ee. 5 

3. Which of the xo tg (i 1,2,...,r-l) in the rth equation is 


to be regarded as Y° and solved for in terms of the other x ts). 


We will call these different ways by the name ''sub-procedures."' 


C. The relation between a remainder and the frequency determinant. 


From the equation 2-4-10, the remainder is 


33 A A 
R(W ) - 33 4&3 + —_— Us, &| 


From equation 2-4-8, we see that A's are the minors of an element of 
the frequency determinant 

Ay = M(U3i) bo= M(Us2) | As = M( Uss) 
where M(u 5, ) means the minor of uz, . 


The remainder is 


A(W) 


es Uz3 A| 


where A(\) ) is seen to be the frequency determinant. 
Expressing a remainder in general in terms of the frequency deter- 


minant; 


R(W ) = Al) (2-4-12) 


Upp coefficient element corresponding to Y° of rth row. 
Arg : minor of u,, in A (LD). 
Upc * element in rth row and the column corresponding to state 


quantity of 1 in Ze: 
22 








D. Failure of the P-M method. 

In case the denominator of the remainder is zero, the P-M method 
certainly fails.* (Refer to equation 2-4-10 and 2-4-12). When the deter- 
minant in the denominator becomes nearly equal to zero or changes sign 
passing through zero, the remainder is difficult to evaluate: at certain 
frequencies the remainder may show an infinite or jumping phenomenon. 

However, this does not necessarily show up at the same frequency 
with different sub-procedures, because this phenomenon depends solely 
on the nature of the transfer matrices. 

Figure 2-4-1 and Figure 2-4-2 show that while one sub-procedure 
shows an inifinite phenomenon at certain frequency, another sub-pro- 
cedure does not show the same phenomenon at this frequency. 

Pestel and Mahrenholtz illustrated their method with the case of a 
homogeneous beam (2r = 4). It leads to very good results for the bend- 
ing frequencies. The reason is, because r = 2, there is no phenomenon 
of indefiniteness causing failure mentioned above and due to purging 
effect. The P-M method does show improvement also for the case of r great- 
er than two. 

For sufficiently high frequencies, the P-M method fails. The reason 
is as follows, if the columns have already become parallel, purging does 
not do any good, and if round-off error reaches critical point, where 


remainder is no longer reliable, this method fails. 


*When this is the case, program VIPIPE skips to the next sub-procedure. 
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Figure 2-4-1 shows fourth and fifth mode frequencies of System 5* at 
frequency 826.88281950 rad/sec, and 1126.83799669 rad/sec (P-M method- 
subprocedure A with double precision arithmetic). Figure 2-4-2 shows 
infinite phenomenon at frequency 939.00298809 rad/sec while Figure 2-4-1 
does not show this phenomenon at this Frequency. And Figure 2-4-2 has 
print out-put with mode frequency 826.88281915, also with 939.00298809 
rad/sec (P-M method-subprocedure B with double precision arithmetic). 
From Figure 2-4-2, we see that the latter is a false mode frequency and 
the former does not show upon the graph due to the infinite phenomenon 
of the latter, but we know that there is a mode frequency at 826.88281915 
rad/sec (this also can be checked with the frequency and remainder value 


print output). 


2-5 Modified Delta method. 

Instead of purging just one column vector among r columns (as in the 
P-M method) of frequency determinant, Marguerre and Uhrig showed in their 
paper [2] a possible way of extending the purging concept for problems 
greater than 4th order. 

Through exactly the same method as Delta, we get a frequency deter- 
minant. At high frequency, because the columns approach parallelism, 
the value of the preoncyadecerittear tends toward zero. If the columns 
are exactly parallel to each other, then of course all methods fail. If 
they become nearly parallel, then with a finite eanmel ty of handling 
numbers, the results become meaningless. | 

Using Maguerre and Uhrig's purging concept, with the help of one 
column of the frequency determinant, obtain r-l purging factors to purge 


*See Appendix E-3. 
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another r-1 columns.* For the case of sixth order (r=3) as used in 
program VIPIPE, the frequency determinant at a certain frequency is: 
A(W)= |YUn Ua Up 
U2 Ugg = Ug 
Uz, Ugg Ugg 


Purging factors got with the help of the third column of A(W) are 


lagigallige eel a) 

, a 3 ( Uy * UB 33 
_ |) eae ee 2-5-1 
dD. = —s( em i! 


Modifying the non-vanishing sub-matrix of the starting state matrix [Zo], 


we have 


0 1 0 (2-5-2) 


Then, for the same frequncy, repeat a round of calculation to get the 
modified frequency determinant. 
Am()= | Uy Ds Up D,Ui3 Us 
Ua D, lbs Ux D, Ub U23 ied,-3) 
Us, D,U33 Us. D, Us3 U33 
From the properties of determinants, A(W) is the same as Am (W). 


But the first and second columns of Am () are not the same as those 


of A(W) any more. If the columns of the original frequency deter- 





minant were exactly parallel, the modified columns (first and second 


columas) of the modified frequency determinant will turn out to be all 


*We use the term 'Modified Delta method'' to refer to this method. 
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zeros. In case the columns were nearly parallel, the first and second 
columns will turn out to be very small numbers, i.e., we will have purged 
two columns. So we can expect a more accurate value of the frequency 
determinant and can go further to higher natural frequencies. 

If D, = D, = O, it degenerates to original frequency determinant, 


1 2 
so it does not do any good. If D, = Dos it likely means that first and 
second columns are the same. If this is the case, this method just fails. 
In some cases, especially systems of many sections, the frequency 
determinant does not show parallelism. Experience shows that this method 


gives good results in most problems. However, this method also fails 


at sufficiently high frequencies. 


2-6 Intermediate conditions - Branched system. 

The effect of a joint point can be introduced into the calculation 
by means of a single point matrix. 

This section is mainly for the derivation of the branch point matrix 
of the Modified Delta method and the Pestel-Mahrenholtz method. 

The derivation of branch point matrix for the straight forward trans- 
fer method (and thus also for the Delta method) is fully explained in 
Ref. [3] and [4] , but a brief discussion will be given here for the 
purpose of refreshing the reader's knowledge and moreover explaining the 
procedure for the M-D method and the P-M method. 

A. Delta method. 

Discontinuities are introduced in the forces whereas the main 
member deflection on each side of joint are the same, and identical to 


those of the branch at the same joint. 
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TO™"SR joint 








= 
0 main member 


branch 


¥ 
Fig. 2-6-1 


By a suitable coordinate transformation, it is possible to express the 
force system of the branch in terms of the main system coordinates. 


A free body diagram for the forces at the joint point is shown 


below.* 
Ma V. Msp Ma. Vv M 
ai JK 
Net ( f |) —7 Nop Na <— C1 \) —> Nor 
Veg Vor 
Va 
v 
oe Nap <— (| 
= / Mae Mna 
Nad 
Fig. 2-2-2A Fig. 2-2-2B 
The forces at the joint in matrix form are. 
M = M an M 
VV V Vi (2-6-1) 


N Ise N Jot N Jas 
The displacement is continuous at the joint, so Z;g, can be obtained by 
simply adding the force components of branch state vector to Z 
The column matrix of nonzero component at the starting boundary will be 


symbolized by V,, Here V, = {1, 1,1} 


*Symbol A indicates that the vector elements are expressed in terms of 
the branch coordinate system. 





~~ Aw 
A 
Ping [Re (2-6-2) 
From equation 2-6-2 
A 7 a 
PeW (RJ(R J + doe (2-6-3) 
The matrix product (RJ- (RJ is called a spring matrix. To transform 


the coordinate system of the branch into main member coordinates at the 
joint, we use 


Ya 


oe = (Game 
7B = (2-6-4) 
Kes = [ G,] - on 
here; G, = (03-0 -sING QO G, “ | 6 () 
SIN D CcOS® 0 0 (0S BG -sin@ 
0 Oo. | 0 swZ tds¢ 
where P is a branch joining angle; refer to Appendix B-5. 
then; 
P ap a} (RaJ- (Ril: (Gi + dx 
= i ’ dag (2-6-5) 


( [ s } = [-&,) ; Gee ac Js (s] is also called a spring 
matrix. ) 


Recalling the relation at the joint 


dsp = 45, = dap 


Pog eee, | Oe We ae Ss [S]- dng 


we get Z,, in matrix form 


he} 


or = [UJ] . 25 (2-6-6) 


Here ([U,] denotes the branch transfer point matrix. 


B. Modified Delta method. 
(In this section, the use of an asterisk * will denote a 
purged matrix or vector. Also note that the quantities d and p, with or 
without * and ~, represent 3x3 matrices whereas in the preceding sub- 
section they represented 3xl vectors.) 
After we get the purging factors from4(W), the new starting matrix 


at location O is 


| 2 
[zJ = j;ol , [N°] = 1 0 
ne 0 l 0 
= = 


At the joint 
] = | a acieg ei) a . (WJ 
P vt 


For a branch, let the non-vanishing state matrix of the new starting 


state matrix [Z,,] be 


[Ng] = 1 @ 0 
0 1 0 
“Dig ~Dag 1 


Here Dj, , D,, are branch purging factors, the calculation of which 


will be explained later. At the joint 


[Z ne *] a4 . (Hh Pe] Le fe 
P dap Ry a 
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[R, ],{R,], (R*,] & [R*,] are defined by this expression. 
Just as in equation 2-6-5, 
(s*) = (@.)-(RI-( R77 (41) 
— Ra}: (Ne) (Nel CRS (61) 
Ga} (Ral (RIT LG 


: os 
The spring matrix turns out to be the same as in the Delta method. 


At the joint 


(CORREIA EV Semi PRA) Ve 
[Pa.J°Vo + [S)>(dnJ-V. 


[Ss] is independent of (No) » so Pe. is independent of [Ne] ~ 


* 
PoR 


* 
For Z sp we have 


Zag® = | Te ©] «28 Dee 
Sa 6 


We see that branch transfer matrix (Up|is the same as in equation 2-6-6. 
Purging factors can be obtained from [R,] or(R}, The phenomenon of 
columns approaching parallel could occur in (R,] and [R,} ; VIPIPE com- 
putes purging factors from [R,]. because the writer was inclined to 
think that getting purging factors from (R,] would do more good in in- 
version of [Rj]. However, [R,J could have been used. 

Purging factors Dip and Do, are obtained in the same way as in 
equation 2-5-1. 

C. Pestel-Mahrenholtz method. 

At the joint the equilibrium condition must be satisfied, 


and the deflections should be continuous. 
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Pag = a Bs Bis 
dag = dn = dep 


At location 0B (boundary of 
a ; 


Here, [Now] is the nonvanishing state matrix of [ 


branch correction vector, 
let 
fo) 0 
[Nc] = P op 0 
Pie | 
1) 
en ¢) 


Then, the state vector at locatior 0B 


is a 
Lop = 0 
0 
0 
Pa Me 
se NE 
Pop Xop X28 
At the joint 
= ; IND XE Ea [47] +X’ 
ae = 
Znp = [4] * (Neo X25 = [Ri (NS) X29= 
| | | 


° 
Nee 


a branch). 


o 





; X oB 
Z°,) and 
0 | Kon = Xe 
Xie 





for the built-in, in-plane 





Xp £8 a 


case 





(2-6-7) 


where [R,*) and[R 


1 *] are 


(RiJ= (Ril (N53) (Ro) = [Ra}'(Noe] 


For(N°} and([X’] , refer to section 2-4. Also [Ry] » [R,] are expressed in 


2 


the same way as the preceding sub-section B. 


(s*] = (@,)-(RE}- (RAJ (61) 

[Gs J+ (Re) (Ne): NF (Rd (40 
(G,J-(Re]-(Ri }-(@1] 

- (3) 


The spring matrix also turns out to be the same as in the Delta method. 


From the requirement for continuous deflections at the joint we get 


= de a 
aye = [6] [RD + (NA) Xs, 
= [ag (0 9. x° (2-6-8 
and 
px = [S)-dae = (8) -[dyf-(n?) - x° 
b>, = ppt x? + [s ]- (ah jx” 
For i we have 
Ze = er ae Taz et TS) ze 
br = 


We see that the branch transfer matrix [Up] is also the same as in 
equation 2-6-6. From equation 2-6-8, the branch correction vector is 
ce) ie +” 4 + . 0 
Xp 7 CR, J06,)-([45]-% (2-6-9) 
After we get the main member correction vector (x ° ), we can get branch 


state correction factors in terms of X° values. As explained in section 


2-4, iteration is performed in the same way. 


33 


2-7 Computer. 

In FORTRAN 63, floating point quantities have an exponent and frac- 
tional part. In single precision arithmetic, the fractional part has 
36 bits and the precision is approximately 10 decimal digits. In double 
precision arithmetic, the fractional part has 85 bits and the precision 
is approximately 25 decimal digits. In both cases the range of numbers 
1#vfven 107°" tou0??. 

In the FORTRAN 63 version of VIPIPE (as described later), either 
single precision (for conservation of computer time) or double precision 
(for greater computational accuracy) can be selected. In the FORTRAN 60 
version, only single precision is available. 

The natural frequencies of a planar piping system are found using 
program VIPIPE. Program VIPIPE has been tested with FORTRAN 60 and 63 
using digital computer CDC 1604. The number of natural frequencies that 
may be found using VIPIPE and the accuracy, especially at high frequency, 
depends upon the floating point significant figure capacity. And the 
significant figure capacity required for any system is a function of 


the characteristics of the system itself. 
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CHAPTER IIT 


DISCUSSION 


3-1 Accuracy of methods and assurance of the solution. 

Confidence in accuracy of methods was established by comparing 
frequencies of several systems for which these frequencies could be 
obtained by other means, i.e., by means of either a closed form analyti- 
cal solution or by recourse to the Ref. [5] . 

The accuracy is very good for the lower natural frequencies in the 
values got using the Delta method, the Modified Delta method, and the 
Pestel-Mahrenholtz's method. 

In most cases, up to fourth mode frequencies, the difference be- 
tween the system natural frequencies obtained by analytical means and 
those obtained using VIPIPE are less than 0.004%. As frequency in- 
creases we got more accurate natural frequency values from the M-D 
method and the P-M method comparing to the original transfer method 
(Delta method). 

In the sample problems reported in Appendix E, we have found that 
for straight, simple configurations (where the phenomenon of the columns 
becoming parallel is pronounced) the P-M and M-D methods give one or two 
more natural frequencies than can be obtained by the Delta method or 
the original transfer method. However, in more complicated configura- 
tions it seems to be possible to get three or four more natural fre- 
quencies. In the simpler cases, there are comparison ail nee available 
from "exact" theory so that the accuracy of the results may be establish- 
ed. However, in the more complicated cases, such comparison values are 
not avelveb1é. However, we have reasonable assurance of the integrity 
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of these results since we have been able to get the same values from 
the M-D method and from several different P-M subprocedures. Thus, 
using these methods, we seem to be able to double or triple the fre- 


quency range in which valid results may be obtained. 





Details of the systems analyzed and the results obtained may be 


found in Appendix E. 


3-2 Limitations and some suggestions. | 
Non-stiff hangers may be introduced into the system by a point 
matrix with appreciable linear and/or torsional spring compliances. 
But, for stiff hangers the transfer method fails, so we cannot use 
program VIPIPE for natural frequencies of a system which has stiff 
hangers. In particular, the suggestion of Fink (4} that rigid hangers 
can be successfully handled by use of idealized very stiff exterior 


springs proves to be untenable. 





Marguerre and Uhrig suggested possible way to overcome the difficulty 
relating to stiff Hanede by introducing unknowns [2] . Introducing un- | 
— means method abandoning the transfer concept, 80 we cannot use 
VIPIPE directly to solve such a problem. However, the writer thinks 
that there is a possibility of modifying VIPIPE to solve such problems | 


by getting state quantities at each node and assuming values for the 


unknowns. For details refer to reference [2] . 
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APPENDIX A 


DESCRIPTION OF PROGRAM 


A-1 General remarks 

Program VIPIPE is a CDC FORTRAN 63 language digital computer program 
designed for use in the vibration analysis of planar piping systems, 
originally developed by George E. Fink (4) . The program VIPIPE is 
modified and augmented by the writer to use Delta method (same results 
as original transfer method, but with less amount of calculation), Modi- 
fied Delta method, and Pestel and aiheamaenel ‘ remainder method with 
double precision or single precision arithmetic. 

In- and out-of-plane undamped natural frequencies may be found 
through the use of program VIPIPE. How this program may be used is 
fully explained in Appendix B. Considerable flexibility with regard to 
the mathematical model employed is provided in the program. For in- 
stance, one may use combination of a distributed mass treatment for 
straight section with a lumped mass treatment for curved sections (this 
method gives the greatest accuracy coupled with least machine time) or 
a lumped mass for all sections. Also one may consider or neglect the 
effect of shear deflection and rotational inertia, singly or together. 

Codes for the boundary conditions and various methods are inter- 
preted within the program. 

Problem solutions proceed by iterative processes with radian fre- 
quency the independent variable. The sign of the frequency determinant 
or remainder is used to control the search for and convergence to a 


system natural frequency. 


Output and transfer to the next method/or transfer to the next pro- 
blem or to the end is provided in the program when it reaches signifi- 
cant figure limitation or previously set iteration limit, and also when 
a method fails. 

All subroutines are provided in binary deck in single precision 
and double precision arithmetic in the FORTRAN 63 version. Modifications 
of the program may be accomplished by . changing FORTRAN cards 
only in the main program; the subroutines are in binary to conserve 


space and time. 


A-2 Program structure 
The program consists of a six section main body and twenty-seven 
subroutines. In the FORTRAN 63 library, GRAPH subroutine is not avail- 
able. So GRAPH subroutine is provided in binary, and built into the 
program. Also, in the FORTRAN 63 library, hyperbolic functions are not 
available, so these functions are computed in exponential form. 
The function of each section of the main body and of each subroutine 
is given below. 
A. Main body sections and their functions. 
INPUT Controls read-in of data, allocates storage loca- 
tions for arrays and subscripted variables, de- 
cides the first method to perform. 

INVARIANTS Computes and assigns appropriate subscripts to such 
system component characteristics as moment of in- 
ertia and radius of gyration, computes frequency 
increment for the first two modes, and sets the solu- 


tion acceptability criterion. 
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CONTROL 
IN-PLANE 


CONTROL 
OUT -OF - PLANE 


ITERATION 


OUTPUT 


Constructs transfer matrix of the system for the in- 
plane case by calling for various subroutines in a 


sequence governed by the system geometry, the mathe- 


matical model specified, and the desired assumptions. 


Performs the same function as CONTROL IN-PLANE but 
for the out-of-plane case. 

Evaluates the frequency determinant (modified fre- 
quency determinant in the M-D method) or the re- 
mainder and utilizesits sign in controlling the 
iterative search for mode frequencies. Initiates 
output when the specified frequency range has been 
fully explored, the required number of modes found, 
the specified iteration number has been reached, or 
the significant figure limit of computer reached. 
Controls output format, provides graph output, seeks 


the method to perform. 


B. Subroutines and their functions. 


SUBSEC 


MATMUL 


FINMAT 


Subdivides components to be treated by lumping mass 
on the basis of L/D ratio and starting frequency or 
required mode numbers. (L/D ratio for curved section 
means ratio of radius of curvature to diameter and 
included angle of arc, as applicable.) 

Constructs the transfer matrix of those components 
treated by lumping mass for only one point mass and 


massless subsection. 


Constructs the system state matrix (6x3) by successive 


premultiplication of the partial system state matrix 
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by the transfer matrix of each component as soon as 
each of the latter is constructed. For the lumped 
mass system, gives successive premultiplication by 
the transfer matrix constructed by MATMUL as many 
times as the number of lumped masses. 

FINBRA Constructs the state matrix of a branch system by 
successive pre-multiplication of the partial branch 
state matrix by the transfer matrix of each component 
of a branch as soon as each of the latter is eenerreee- 
ed. For the lumped mass system, gives successive pre- 
multiplication by the transfer matrix constructed by 
MATMUL as many times as the number of lumped masses. 

STATEM Constructs starting state matrix from the starting 
boundary condition of main member, for each method. 

STATEB Constructs starting state matrix of a branch system 
from the far end boundary condition of a branch for 
each method. Computes the purging factors of the 
branch in the M-D method. 

BRCOR Constructs branch correction matrix for the assumed 
branch starting state matrix at a certain frequency 
in the P-M method. Constructs branch purging matrix 
in the M-D method. 

DELMA Constructs frequency determinant in the Delta method. 
Evaluates the purging factors from the frequency 
determinant of the Delta method and constructs the 
modified frequency determinant after another round 
of calculation in the M-D method. For the P-M 
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method, computes the correction factors for each sub- 
procedure. For each sub-procedure, if it fails by 
the fact that denominator is zero, the rows are 
changed and a second test is made. It also provides 
transferring to the next method or sub-procedure if 
it fails. 
DISTM, DISTMO, SFIELD, SFIELO, CFIELD, CFIELO, POINT, POINO, RIGID, RIGIO, 
STIFCO, HANGER, HANGEO, BRANCH, BRANCO, STAVEC, STAVEO, INVERT. 


For the above subroutines refer to reference [4] . 


A-3 Program nomenclature. 

B1,B2,B3 Correction factors for assumed state quantities of 
the far end of a branch. 

BC Discriminant. Controls boundary correction in the P-M 
and M-D method. 

BL Number of the last iteration. 

D1,D2 Purging factors in the M-D method. Correction factors 
in the P-M method. 

DEL Array of numerical values of frequency determinants 


and array of remainders. . 


DV Numerical value of a determinant in denominator of 
X1 and X2. 
EDL Array of numerical values of frequency determinants 


or remainders in single precision for input to GRAPH. 
ERM A discriminant same as REM in the problem statement 
data card. 


FF Discriminant. Controls output to provide a print 
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output that analysis is terminanted due to signifi- 
cant figure limitation in the P-M method. 

FINK Discriminant. Controls constructing the frequency 
deeerninane or modified frequency determinant in the 
M-D method. 

HH Discriminant. Controls transferring to another sub- 
procedure or to another way of testing in case one 


fails in the P-M method. 


HMH Discriminant. Controls subsectioning in lumped mass 
7 system. 

K | Discriminant. Controls grid in the GRAPH output. 

KA Discriminant. Terminates iteration when specified 


iteration number is finished. 

KKK Discriminant. Controls transferring to the output 
in case the last sub-procedure in the P-M method 
fails or controls transferring to another three sub- 
procedures (D,E,F) after three sub-procedures (A,B,C) 
are finished in case subprocedure C fails in Selec- 

‘.tion 3 (for Selection, see Appendix B-9-6). 
MAT Discriminant. Controls matrix multiplication in the 


lumped mass system. 


NUMPTS Number of points in one graph (except the last). 
NUMEND Number of points in the last graph. 

res PZ Assumed starting state quantities of the main member. 
PMM Discriminant. Controls three methods (Delta, M-D, 


P-M method) specified. 
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PPM Discriminant. Assures that same procedures are used 


for out-of-plane solution that were used for in-plane 


solution. 
PM Discriminant. Controls three methods to be used. 
PP Discriminant. Controls the three sub-procedures 


(A,B,C) in the P-M method. 

RR Discriminant. Controls testing by interchanging two 
rows in case one sub-procedure fails because the de- 
nominator of Xl or/and X2 is zero in the P-M method. 

R3 Array name. Matrix product (r,) “ ([¢,] , used in 
the branch correction. 

REM Discriminant. Controls three ways of selecting sub- 
procedures in the P-M method. 

Q1,Q2,Q3 Assumed starting state quantities of a branch. 

ss Discriminant. Controls the three sub-procedures (D, 
E,F) in the P-M method. 

XX Array name. Branch correction matrix (3x3) in the 
P-M method and branch purging matrix in the M-D 
method. 

X1,X2,Y¥1,Y2 Correction factors for assumed starting state quanti- 
ties of the main member in the P-M method. 

YY Array name. Branch correction matrix (3x3x23). 


{ 
For other nomenclature’ see reference [4] : 
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APPENDIX B 


INSTRUCTION FOR PROGRAM USE 


B-1 General remarks. 

Appendix B is intended to provide complete instructions for using 
program VIPIPE. VIPIPE is FORTRAN 63 language program with twenty-seven 
subroutines. All subroutines are provided in a binary deck so that the 
compiling time is reduced as is the volume of the deck. 

Program VIPIPE can be see single precision or double precision 
arithmetic. Also VIPIPE is provided in FORTRAN 60 language with single 
precision. 

The Delta method, the M-D method, and the P-M method are all incor- 
porated in program VIPIPE. The M-D method and the P-M method do appear 
to permit obtaining more natural frequencies (one,two, or possibly several 
more) than can be obtained with the Delta method (same result as the 
original transfer method). However, the P-M method has the disadvantage 
of introducing an odd behavior of the remainder. In some cases, this 
odd behavior can easily be misinterpreted so as to lead to false fre- 
quency determination. However, there are a whole group of associated 
P-M modifications* and the odd behaviors occur at different frequencies 
for the different modifications. Thus, by making the calculation more 
than once, using a different modification each time, the true natural 
frequencies can be sorted out. 

Program VIPIPE uses the Delta method as the basic procedure. If 


more frequencies are required than have been obtained by this method, 


*P-M modifications means sub-procedures of the P-M method. 
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one may then use one or more of the available associated P-M modifica- 
tions and/or the M-D method. Graphical output will assist in identifying 
the natural frequencies. The writer recommends to use the P-M modifica- 
tions as original method in getting more natural frequencies, and to 

use the M-D method for purposes of comparison. 

The accuracy and number of natural frequencies that can.be obtained 

using the methods depends on the capacity:of handling numbers accurately. 

As we can see in Appendix E, double precision arithmetic gives better 

results especially for complicated systems for which many matrix multi- 
plication are performed. But using double precision arithmetic results 

in considerable increase of machine time. And also for various methods, 

the machine time is different. The Delta method and the P-M method each 

take about the same amount of machine time, but the M-D method takes : 
almost twice this time. 

So a user should consider the frequency range under investigation and 
the nature or characteristics of the system as well as the machine time. 

For comparison purposes Section B-9-6 of this Appendix lists machine 
times for one iteration for some example problems using each of three 
methods with double and single precision arithmetic. 

Complicated systems are divided into one main member and branches. 
System components are defined as sections of pipe, straight or of con- 
stant curvature, an elbow, a spring hanger, a coupling, a flange. Dia- 
meters, wall thickness, density, elastic and shear modulus may vary from 
component to component, but may not vary within a component. Both tors- 
ional and linear spring rates (for hangers) must have finite values. 


Thus perfectly rigid hangers cannot be treated. 
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B-2 Coordinate system. 

A local coordinate system is constructed at each point of interest 
in the system. A right handed coordinate system is used whose X axis 
is the local tangent to the central axis of the main member, always 
directed away from the left end. The X-Y plane coincides with the 
quiescent plane of system. The Z axis is always directed downward. The 
direction of the local Y axis is chosen so as to complete the right hand- 


ed system. 





Fig. B-2-1 


The branch coordinate starts from the end of the branch farthest from 
the joint point with the main member. The plane of the branch coincides 
with the X-Y plane. The positive Z axis of the branch and main member 
are parallel and have the same sense. The Y axis is oriented so as to 


complete a right handed system. 


B-3 Data required. 
Extensive properties. 
a. outside diameter (D) - inches. 
b. wall thickness (DI) - inches 
c. radius of curvature (RHO) - inches (Negative, if center of 


curvature lies on negative Y axis). 
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j. 
k. 


Intensive 


B-4 Component 


length (TL) - inches (for straight section only; for 

curved sections, length is computed by the program). 

arc central angle (THETA) - degrees 

linear spring constant in X direction (CLX) - pounds per 
inch 

linear spring constant in Y direction (CLY) - pounds per 
inch 

linear spring constant in Z direction (CLZ) - pounds per 
inch 

torsional spring constant about X axis (CTX) - in-1lb/radian 
torsional spring constant about Y axis (CTY) - in-1lb/radian 
torsional spring constant about Z sxis (CTZ) - in-1lb/radian 
properties 

density of component (AAMU) - lb/ft? 

elastic modulus (AE) - psi 


shear modulus (AG) - psi 


identification.* 


Components are identified by two numbers. 


A. All components have a component sequence number starting from 


the most left component of the main member. This sequence number is 


identified by the order of data cards containing component extensive and 


intensive properties. This can be understood from the figure B-4-1. 


*Aslo see pp. B-5 to B-7 of reference [4] . 
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Fig. B-4-1 


B. A branch component identification number (NNBR) is read in 
as data. It indicates whether it is in a branch or main member, if it is 
in branch, the position in branch. 
NNBR O main member 
NNBR 1 the first component of a branch (farthest from the 
intersection) 
NNBR 2 intermediate component of a branch 
NNBR 3 the last component of a branch 
NNBR -1 the only component of a branch 


This can be understood from the figure B-4-2 






aoa Y oi y 


S 2 
Branch 3 





Branch | 


Fig. B-4-2 
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B-5 Branch joining angle 

Branch joining angle to main member is measured counter-clockwise 
in degrees from the X axis of the main member, as it is oriented immedi- 
ately before the inter-section, to the X axis of the branch as it is 


oriented at the point of intersection. 


B-6 Boundary conditions* 
Five common boundary conditions, those of a fixed, free, pinned, 
propped, or a roller supported end are indicated by the code given. 
The program interprets the code and forms the required state 
matrix (6x3). 


Boundary condition code; 


a. sv l fixed end 

b. SV 2 free end 

c, SV 3 pinned end 

d. SV 4 propped end 

e. SV 5 roller supported end 

f. SvVoO other than one of the above. 


(The user must construct the necessary state 
vector and cause to be read in as data). 


This code is used for both in-plane and out-of-plane cases. 


B-7 Control cards. 

A user sets up for execution with a master control, a FORTRAN con- 
trol card, and combination of END, FINIS, EXECUTE, transfer cards, bin- 
ary cards, and FORTRAN cards properly placed in FORTRAN 63 deck (JOB and 


END cards are used in FORTRAN 60). Figure B-7-1 shows the deck assembly 


*For further details, see pp. B-8 to B-10 of reference [4]. 
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for FORTRAN 63 use. 








he) 


é 
v 


= a ” 
. 


DATA 
: DECK 
A 
We BINARY DECK 
SUBROUTINE) 






PROGRAM VIPIPE 





§ COOP, ,KIM Y S,0/49/S/1S/2S/E/45 54,30,30000. 
FORTRAN DECK 


(MAIN PROGRAM) 


Fig. B-/7-1 


B-8 Data deck assembly. 

The data deck is ordered as follows: 
a. problem statement card. 
b. boundary condition data card(s). 
ca branch joining angle data card(s). (deleted when system 

has no branches) 

ly component extensive properties data card(s). 
e. component intensive properties data card(s). 

The user may obtain solutions for any number of problems desired; 


in this case each problem data deck is assembled just behind the other. 
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Following the last problem deck a blank card should be added. 


B-9 Flexibility of program VIPIPE,. 
B-9-1 Double precision and single precision arithmetic. 

Program VIPIPE can use double precision process only with FORTRAN 63. 
Double precision is used in system state matrix (UU), branch system state 
matrix (VV), and frequency determinant or remainder. 

Type double is declared in main program and subroutines FINMAT, FIN- 
BRA, STATEM, STATEB, BRANCH, BRANCO, BRCOR, DELMA. By removing cards 
declaring type double the program can easily be changed into single pre- 
cision. Subroutines are provided in binary deck with single precision 
and double precision arithmetic, so with corresponding binary subroutine 
deck and main program in FORTRAN with or without a card (card number 
000120) declaring type double we can use single or double precision arith- 


metic. 


B-9-2 GRAPH. 

A user can get GRAPH output by option. The values of the frequency 
determinant and the values of the remainder are not usually of the same 
order of magnitude throughout the frequency range. Moreover, because 
of the infinite phenomenon of the P-M method, it is not generally useful 
or convenient to get one graph for the whole frequency range. 

1. Number of points in one graph can be changed by changing 
NUMPTS in card 006090. 
2. Can get graph with grid or without grid by changing card 


006100. 


aw 
il 


1 with grid 
K = Q without grid. 
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3) Graphs are labeled from 1,2,3,...up to 10. Graphs beyond 
tenth one Will pe labeled as 10. The last graph is label- 
ed LAST. 

4. For each graph, the whole frequency range is moved (except 
the first one) so that the first point of each graph starts 
from 0. The corresponding value of frequency determinant 
(or remainder) for the first point and NUMPTS (NUMEND for 


the last graph) are printed in the print output. 


B-9-3 Shear area form factor. 
A shear area form factor of one half is built in program.* This may 


be changed by changing card 000600. 


B-9-4 Starting frequency, frequency increment and solution acceptability 
criterion. 

Using the M-D or the P-M method, we are expecting to get more natural 
frequencies at high frequency. Thus choosing the proper starting fre- 
quency for the M-D or the P-M method, following the Delta method can 
save machine time. This can be done by putting the wanted starting fre- 
quency in problem statement data card. 

To compute the starting frequency increment, the program constructs 
a synthetic straight pipe, equal in length to the length of the main 
member, and having diameter, wall thickness, and intensive properties 
equal in magnitude to a weighted average of those properties of the 
composite members. One quarter of the fundamental frequency of this 
synthetic pipe is taken as the starting frequency increment (card 001390). 
And one thousandth of this increment is taken acceptability criterion 


(card 001410). When the iteration reaches this solution acceptability 


*See reference [7] . | 53 


criterion, the corresponding natural frequency is determined using linear 
interpolation. 

This frequency increment is for the first two mode frequencies, but 
the acceptability criterion is the same for all mode frequencies. After 
second mode frequency is found, frequency increment is one tenth of the 
difference between the last two mode frequencies found (card 003950). 

These can be altered by changing cards. The writer found that a 
smaller frequency increment gives better results especially at high fre- 


quency with the P-M method. 


B-9-5 Subsectioning. 

Subsectioning for the lumped mass treatment of components is carried 
on in subroutine SUBSEC. A component is divided into maximum of twelve 
subsections. 

For all components having a length to diameter ratio greater than 
six, the number of subsections is computed as two times the value of a 
parameter HM which appears in SUBSEC (HM is number of modes to be sought). 
If HM is greater than six, HM is given the value of six. In case the 
starting frequency is greater than 800 radians per second, HM is also 
given the value of six. HM can be controlled by a parameter HMH in the 
main program (card 001250). If HMH is zero, HMH does not affect the 
value of HM. By changing HMH to other than zero we can change the value 
of HM equal to that of HMH. 

Poor subsectioning might give some loss of accuracy in determining 
natural frequencies, especially at high frequency. Checking the natural 
frequencies by increasing the number of subsections using parameter HMH 


is recommended. 





B-9-6 Machine time and printing line limit. 

In FORTRAN 63, estimated time is set to thirty minutes and print- 
ing line limit is set at 30000 lines in VIPIPE. These limits can be 
altered by changing the master control card. If a limit is reached, 
computation stops automatically. 


For reference, machine time for one iteration of some example problems 
is given below. 














Double 
Prec. 



















B-9-7 Performing various methods. 

One selects the methods to be employed by use of the problem state- 
ment data card. Then the program picks the first method in the INPUT 
section of the main program with a discriminant PM, in the order: Delta 
method, M-D method, P-M method. 

Finishing a method, the OUTPUT section of the main program opens to 
the next method. Finishing all specified Mineas, the program goes to the 
out-of-plane portion, to the next problem, or to the end. 

In the P-M method, there are six sub-procedures. These sub-pro- 
cedures are controlled by discriminatns PP, SS. 

PP.: 0. subprocedure A, SS : 0. sub-procedure D, 
Beaee lL. sub-procedure B, SSac: lk. sub-procedure E, 
Pre: 2. sub-procedure C, SS. 3 We sub-procedure F. 


*For example systems refer to Appendix E. 
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Each sub-procedure has another possible way by interchanging rows of final 
determinant (refer to section 2-4). If one fails, the alternative one 
is tested. If one succeeds, then go to the next sub-procedure. There 
are three ways of selecting sub-procedures by discriminant REM. 
REM: 2. three sub-procedures in the order: A,B,C. - Selection 1 
REM: 1. three sub-procedures in the order: D,E,F. - Selection 2 
REM: 0. six sub-procedures in the order: A,B,C,D,E,F. - Selection 3 
If one wants just one sub-procedure, by changing PP to a value of 2 
(card 000710) in selection 1 only sub-procedure C is used, or changing 
SS to a value of 2 (card 001000) in selection 2 only sub-procedure F is 
used. If one wants just two sub-procedures, by changing PP or SS to a 
value of 1 sub-procedure B and C or sub-procedure E and F, respectively, 


are used, 


B-9-8 Limiting the number of iterations. 

One can limit the number of iterations, for each problem, by a para- 
meter HL in the problem statement data card. In the P-M method, because 
of the infinite phenomenon, the frequency increment in iteration can be- 
come very small. Thus the iteration might keep going on indefinitely. 

The program permits six hundred iterations as maximum. The writer 
found three hundred iterations to be adequate to get nine natural fre- 


quencies with the acceptability criterion set in the program. 


B-9-9 End conditions. 

A boundary condition state vector must have three non-zero and three 
zero elements. Frequently, however, an actual condition may not meet 
this specification. For instance, an end may be supported or connected 
to a piece of machinery. Then we approximate it as being held in a 
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mounting block having specified linear and torsional spring rates, and 
to this attach a fictious massless short section. 

Through the use of this sort of fictious end, a proper end condi- 
tion state vector may be constructed for any real situation with little 


or no loss of accuracy. 


For hangers exhibiting coupling, systems having unusual flexibility, 


and flanges and couplings refer to pp. B-19 to B-20 of reference (4) . 


B-10 Format of data 
B-10-1 Problem statement 
The card has 15 fields under control of the following field specifi- 
cations. 
312, 4F3.0, 2F8.0, 4F2.0, 2F4.0 
Field names and the information each field conveys are given below in 
order. 


A. Field names and content. 


1. NS number of components. 

2. #$NBR number of branches. 

3. IOP code for in-and/or out-of-plane solutions. 

4. HM number of modes sought. 

5.  AMYK code calling for the use of lumped or distri- 
buted mass model 

6. SD code calling for shear deflection inclusion or 
neglect 

i RI code calling for rotational inertia effect in- 


clusion or neglect 
8. OMEGAG upper limit of frequency range of investigation 
9. OMEGAT lower limit of frequency range of investigation 
and initial iteration frequency. If 0 is placed, 
it will begin at .01 radians per sec. 
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10. 


Lig 


LD. 


16. 


GDABC 


SD 0. 


RI 0. 


RI 1. 


GDABC 0. 
GDABC 1. 
GDABC 2. 


BRC 0. 


BRC l. 


GRAPH 0. 


GRAPH 1. 


code calling for print output... 

code calling for the inclusion or exclusion of 
branch correction factors in P-M method and 
purging factors in M-D method. 


code specifying whether there is to be graph 
output. 


code for choosing methods in P-M method. 


code calling for the combination of 3 methods of 
Delta, M-D, P-M method which a user wants. 


number of iterations. 


Both in-and out-of-plane mode frequency sought 
(in this case, in-plane frequencies are sought 
first). 

In-plane frequencies are sought. 

Out-of-plane frequencies sought. 

Lumped mass model, 

Distributed mass model, 

Consider shear deflection. 

Neglect shear deflection. 

Consider rotational inertia. 

Neglect rotational inertia. 

Print solution only. 

Print solution and input data. 


Print solution and input data and graph data. 


Do not consider branch correction factors or 
purging factors. 


Consider the branch correction factors or purg- 
ing factors. 


No graph output. 
Graph output. 
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17. REM QO. All six sub-procedures in the P-M method (A,B, 
CU shar 


18. REM 1. Three sub-procedures in the P-M method (D,E,F). 
19. REM 2. Three sub-procedures in the P-M method (A,B,C). 
20. PMM l. Delta method only. 

21. PMM 10. M-D method only. 

22. PMM 100. P-M method only. 

23. PMM ll. Delta, M-D method. 

24. PMM 101. Delta, P-M method. 

25. PMM 110. M-D, P-M method. 


26. PMM 111. Delta, M-D, P-M method. 


B-10-2 Boundary conditions. 

Field specification 22F/.3 

Boundary conditions are read in separately for in- and out-of-plane 
cases. When both cases are to be investigated, in-plane case data are 
read in first. 

N + 2 fields must be filled in (N is number of branches). The first 
and last fields are first and last boundary condition codes of the main 
member. Branch boundary conditions are entered in the intermediate fields 
in accordance with the branch sequence number. When a boundary condition 
is not describable by one of the codes, a zero must be entered in its 
data field. For each zero an additional data card is required. The 
limitation in non-zero state quantitites at the boundary must be three 


(refer to B-9-9). 


B-10-3 Branch joining angle (degrees). 


Field specification is: 20F7.3 
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B-10-4 Component extensive properties and branch component identifica- 


tion. 
A data card for each system component, ordered by sequence number. 
Field specification is: 3F5.2, 2F6.2, 6F7.2, 12 
Field names are: in order, 


D, DI, RHO, TL, THETA, CLX, CLY, CLZ, CTX, CTY, CTZ, NNBR 


B-10-5 Intensive properties. 


The first intensive properties data card has four fields with speci- 


fication: 


F7.3, 2meo.2, "2.0 


Field names are: AAMU, AE, AG, DD* 


*When the intensive properties of all components are the same, the DD 
field is left blank. When the intensive properties are different, the 
DD field is non-zero, then an intensive data card must be prepared 
for each of the remaining components (including hangers). Additional 
data cards are identical to the first one except that the DD field is 
blank. 


A0 





APPENDIX C 


FLOW DIAGRAMS 


C-1 General remarks. 

Flow diagrams are included for all parts of the main program, and 
for several subroutines. Since CONTROL OUT OF PLANE, FINBRA flow dia- 
grams are identical in structure to those of CONTROL IN PLANE and FINMAT 
respectively, only the latter are included. Reader can refer to the 
flow diagram of STAVEC and STAVEO of reference [4]. Since the flow 
diagram structures of SUBSEC, BRANCH and BRANCO are similar to those 
of reference [4], they are not included here. 

Flow diagrams of other subroutines that are not included in this 
paper are relatively simple ones, so the reader can refer directly to 


the program. 


C-2 Index of flow diagram. 


Main program. 


Flow Diagram. Page. 
INPUT 63 
INVARIANTS 65 
CONTROL IN PLANE 66 
ITERATION 68 
OUTPUT 70 
Subroutines 
FINMAT TZ” 
STATEM 74. 
STATEB T5 
BRCOR 76 
DELMA 77 
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C-3 Flow Diagrams 
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INPUT FLOW DIAGRAM _ | — = . 





READ: Dimension 
; statement 
















READ: Problem 
statement 











70 
N 
READ: Coded case2 READ: Coded case 1 
boundary conditions boundary conditions z 


tea 


: r 
READ: Coded state 


vector SVOP(I,J) 
T=1,6 







READ: Coded state 
vector SVIP(I,J) 


READ: Branch joining 
angle PHI(I I=1,NBR 







READ: D(I), DI(I),RHO(I), Tiglela) THETA(I) 
Cinta) CLY(I), CTZ( ae crx(T), cTy(I), 
CYZ(I},NNBR(I) Im1,NS 


READ: AMU(I),E(I),G(I) 
|e T 260 


OMEGAI=.0 


3 ty 
On + 
(yu z SE 
ny 8: 
> 
w 
nS is 


INPUT FLOW DIAGRAM(CONT., ) 


AMYK = AMYK-1. 

OMEG AO = OMEGAT 

FFAC =.5 . 

L=1 

M=1 

AA™=BB=CC =GG=EE=0, 
AAMU== AL== AD =AID=0O. 


a 


KK=KKK=RR= PP =KA=0 











>0 
=f <Htr0 01> () 
PM = 1, 


THETR(I) = THETA(T) 
DT(I)= ih 

BAMU(I) = AMU(I) 
DI(I) = DI(I)-2.0*DI(I) 
I#x1,NS 






GO TO INVARIANTS 


Li Ww G 


FROM INPUT 









CLX(I)=CLY(1I)=CLZ(I) 
CTE(1 )msieTY(1)=cTra@) 
INS 


=60 


nae = anu(I) *Ee(D(1) ** 2-DI(1)**2) (69120 *32,17*12. ; 
THETA(I)= THETA(I)*7/180. |#40% 
TL(1) = THETA(I) *ABSF(RHO(I 


Ag (I)= W*(D(1) **4D1 (I) **4)/32. 
AJ(1)=aJT(1)/2. 

R(T )=tl Case ner) 

AIX(1) = SQRTF((D(1)**2 +DI(1)**2)/8 
ALY(I) = AIX(1I)/1.414214 


- 
| 





HHM—HM BMH=0. 
CALL SUBSEC 





AD = AD+ D(1)*TL(I) 

EE= FE+E(1I) *?L3(1) 

AID= AID+ DI(I) *TL(I) 
AAMUs= AAMU f- AMU(1)*TL(I) 













AD = AD/AL 

AID = AID/AL 

AFR =1764.*( AD**4-AID**4) | 
ADOM = FR/4.0 

DOM = ADOM 

CR= FR*0.00025 
FR = 4.73** 2*SQRTF( EE*AFR/( AAMU*AL**4 ) 


BDOM = ADOM 
BR= CR 





6a 


ENVARIANTS FLOk DIAGRAM(CONT.)__ 


HI(IT) = PHI(II)*%/180. 
II =1,NBR 






O CONTROL IN PLANE 


TO CONTROL OUT .OF PLANE 


CONTROL IN PLANE FLOW DIAGRAM 
FROM_INVARI AN 
Pl= P2=Q1= Q2=Q3 =1 


D1=f2=Bl =B2.=B3= 0 
FINK= 0. 


SECS ail STP) —— en (CALE STAVEC 
Jm + 2 ms 8, 


SVIP(J,NMBR)}e——* 


si a 
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CONTROL IN PLANE DIAGRAM (CONT. ) 


ee ee ey age tears) er. 


. P( I) = CLX(1)-+ CLY(1)+ CL2(1) MAT= 0 





ALL RIGID; |CALL DISTM 


CALL BRCOR 


=} 
TO ITERATION 


—— ee I 





oe me 


Numericel velue of | 


frequency determinant 


“GG <2, as 
HDEL = DEL(L) 
OMEGA(Irr1) = OMEG2(L)+ 2.*CR 


a ae at Alin oe, 


TERATION FLOW DIAGRAM 





CALL DELMA 


70 ) 


GO TO: Address 
6001 of OUTPUT 





DEL(L)=™ X1-Y1 
Di=(x1t Y1)/2. 


D2=X2 





DEL(L)> 


>0 


x 
~ 


<( 
(oy DELC* DEL(L}> 


70 


=s OMEGA(L41)= OMEGA(L)- 


40 ey 
<Co> {oMBG A(L) -“oMBCAG PDEL = DEL(L)_ 






PDEL DEL(L-1 


. (OMEGA(L-1) -oMBGA(L)) *ABSF( DEL(L 
one: ¥i(M) =CMEGA(L) (oupea( te -owesa(s > ABSP( Pett) 
A(M) = CHEGAC “/¥ sBSR(DEL(L))-+ ABSF(DEL(L-1 


es OS SS Ee Ne mE ee > 2 wae 


(PR 





ITERATION FLOW DIAGRAM( CONT. ) 


DELC = DEL(L) z 0 | 
OMEGAC = rua yt as BB=O. PDEL*DEL(L 


ee We. 0 <3 
@®oneca(t a sO GrBGR L)-OMEGAG)4M=M-1] <A> 


|70 OUTPUT | 


OMEGAC = OMEGA( L-1) 







DELC = DEL(L-1) 


‘ 


ee oe >0 
Joweca(at2) = omeca(L) -(onaea(L)-onmcac) /2. <p> 






£0 /omeca(L)-OMEGA(L-1)-CR 
aes 













C 
= % 
“BCAM(¥)= 0 OMEGAC+ c+ Loma fe) ana =e Orin OMEGA(L)-OMEGA(L-1)-CR 
>0 | 
DELC =0. 





AA= BB =CC =GG=0 






OMEGA(Irve1) = OMEGA(L)- 
DOM= .5*DOM 


- 9* DOM 








GO TO Addressl002 
CONTROL OUT OF PLANE 
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OUTPUT FLOW DIAGRAM _ 


° ——eweeeew © = © © e © 8 FF 


FROM _I'' ERATION 


a ee 0 
PRINT: Heading /{#—<(> 
=<ID 270 


PRINT: IN-PLANE PRINT: OUT-OF-PLANE 
fresuencies frequencies 


PRINT: Mathmatical 
model used 











0 
PRINT: Deltas PRINT: M-D PRINT: P-M method, 
method used method used 6 subprocedure used 
=f DABC - | ? aloes Derren 





PRINT: Subsection data 
lgraph data,input data 


PRINT: Input data| <0 


Sort OMEGA 
aS, 


ee oe we 







, Address 6001 
“1; 0 


MH=M LM], ‘i 


OMEGAM(M)=0. 
Me 


| ORBGACE)= el 
L=1,1M _ 
t 


DELC = OMEGAC=0. 
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OUTPUT FLOV DI AGRAM(CONT. ) 
















= AA=BB=CC =GG=KK= 0. 
ADOM=BDOM 

CR=BR 

OMEGA(1)= OMEGAI 
OMEGAO=OMEGAI 

DOM=ADOM 

KKK=>KA= FF=0O. 

L=M=1 





GO TO: Address 450 INPUT 
















SELECT: Next 
specified method 





Finished all 
specified method 





GO TO: CONTROL OUT 












SUBROUTINE FINMAT FLOW DIAGRAM | 





v(J,K)=0. 













Si Joi 6 K=1,3 
N= Z-1, . 
a V(J,KV(S,K)+U(J,L)*UU(L,K) 
' iS Je1,6 Kel,3 L#1,6 
V(J, \= 0, 
J 1,6 Kegs 
VWI, K)=V(J,K)FU(J,L) *OU(L, x) 
_J=1,6 Kel,5 Url,6 





uu(s, K)=U(J ,K) 
Je1,6 Kel,3 


20 


V(J,K)=0.. 
J=1,6 Ke1,3 





v(J,K) v(J,K) A(J,L)*UU(L,X) 
J<1,6 K=1,3 L=1,6 


ae I 


| RETURN 


55 


yf: 


SUBROUTINE STATFM FLOW DIAGRAM 






- YUCE, K=O. 
J=1,3 Kal, 3; 


<ep>—<S 


>9 


Pl=PltD1l 
} P2™=P 24D 2 


Me@= Mtl 


— | vu(s,M)=1. 


=0 
M-1> UL J = 
re 3}=Pl 


UU(J,3}D2 


ZO 

= 060 Cee—— 

KU-D uU(J,3)=P2 
#0 | 


meee seme ee 


SUBROUTINE STAT2B FLOW DIAGRAM 





w(J,2)=-Dl 
—> 


CALCULATE: 
Purging factors 
Dl & D2 


J-1,3 K=1,3 


> a?) 


[vv (s M)=1. 





A &D 


J= Jt 


RETURN 


CALCULATE: 
Bl,B2,B3. 
For subproc. 








—-—— oe 


— 


 D1l= D2=0. 










<0 >0 


C : A 
=O | CALCULATE: 






For subproc. 
C& F 





CALCULATE: 
Bape, BS. 

For subproc. 
B& ED 







Ql=Q1*Bl 
Q2= Q2*Bl B2 : <a 
Q3=Q3*Bl B3|. | 








SUBROUTINE STATEB » LOW DIAGRAM( CONT. ) 





vv¥(J,1)=1. 










XX(J,K)=VV(J,K) XX(1,K)=VV(1,K uuu(1,K ; ; 
J=1,3 Kel,3 eae ld ee. a2 a 3,8 Re mag 
PS 47) uuU(3,K)=UU(4,K) 
= K™1, 3. 


XX(J,X)=0. 
J=1,3 K=1,3 


XX(J,K)=XX(J,K)+R3(J,L)*UUU(L,K) 
J=l,3 K=1,3 L™1,3 
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SUBROUTINE DELMA FLOW DIAGRAM 


| bax = k2=-Y1= Y2—D1= D2= 0] 








CALCULATE: 
et | <0 Purging 
SI =O <b 6 =" factor D1l,D2 
70 
Gp <p |” 


= 
C2 
KKK =1 D(2,K)=Vv(3,K) 
D(3,K) = V(2,K) 
K=1,3 
RETURN CALCULATE DV| 





eatin 


See’ 


a 


CALCULATE 
Xl & Xe | 


yi 






CALCALATE Y2 CALCULATE Y1 









[RETURN] 
(A) 
ss) = V(3,K) D(J,K)=V(J,K) RR= 1. 
= V( 2) Je2,3 K=1,3 HH =1. 
152 


CALCULATE DV 





SS=1j (P=. 
CALCULATE DV 


RR=0. 
HH=1. 
OV OY) 
0 
RETU RN i 
i CALCULATE 
XT & e 


FF=1. eg 
RETURN CALCULATE Y2 | CALCULATE 


RETURN 


RETURN 
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SUBROUTINE DELMA FLOW DIAGRAM( CONT, ) 





CALCULATE DV] : 





2 gs = 
RETU RN CALCULATE Y2] [CALCULATE Y1 






RETURN 
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APPENDIX D 


PROGRAM LISTING 


D-1 INDEX 

Main Program. 
Section 
INPUT 
INVARIANTS 
CONTROL IN PLANE 
CONTROL OUT OF PLANE 
ITERATION 
OUTPUT 

Subroutines 
SUBSEC 
DISTM 
DISTMO 
RIGID 
RIGIO 
STIFCO 
STIFOO 
STAVEC 
STAVEO 
INVERT 
HANGER 
HANGEO 
CFIELO 


POINO 


Sequence 


000000 
001080 
001510 
002460 
003400 


004430 


008000 
008530 
009100 
009670 
009930 
010190 
010440 
010690 
011000 
011270 
012010 
012190 
012370 


012890 
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Number 
001070 
001500 
002450 
003390 
004420 


007990 


008520 
009090 
009660 
009920 
010180 
010430 
010680 
010990 
011260 
012000 
012180 
012360 
012880 


013100 


Page 
82 
83 
84 
86 
88 


89 


96 
97 
98 
99 
99 
100 
100 
101 
101 
102 
103 
103 
104 


105 


D-2 


Section 


CFIELD 
POINT 

SFIELD 
SFIELO 
MATMUL 
FINBRA 
FINMAT 
BRANCH 
BRANCO 
STATEM 
STATEB 
BRCOR 


DELMA 


Listing. 


Sequence Number 


013110 
013570 
013800 
014160 
014430 
014570 
014860 
015150 
015790 
016480 
017030 
017790 


018110 


81 


013560 
013790 
014150 
014420 
014560 
014850 
015140 
015780 
016470 
017020 
017780 
018100 


019500 


Page 
105 
106 
106 
107 
107 
108 
108 
109 
110 
111 
112 
113 


114 


ANNAAN 


450 
1500 
601 
1501 
1502 


1503 
1504 
F205 


1506 
1507 


1508 
1509 
1510 
Lo23 
Poi) 
PS 12 
ols 


1514 
Pols 


1516 
ey 7 


1518 


ae 
1520 


PROGRAM VIPIPE 
INPUT 


DIMENSION AJT(50) 29AJ(50) 29R (50) sAIX(50) sALY(50) 9PHI (20) »SVI(601)9 
1DT( 50) sD(50) 901 (50) TL(50) 92150) 36150) 5E (50) » AMU( 50) » THETA( 50) » RHO 
2(50) sOMEGAM( 20) »DEL (600) sOMEGA ( 600) 9A(696) 9B (606) 9U( 606) sCLX(50)0C 
ALY (50) 0 CLZ(50) 9 CTX(59) sCTY (50) »CTZ(5C) 9V (693) 9 VV( 603) 0UU( 693) »SVE( 
4691) 9SV(22) sSVIP(6922) »SVO(22) »SVOP (6922) sNNBR( 50) »BAMU(50) » THETR( 
520) 9P (20) sYY (393923) 9R2(393) oXX( 393) eT TITLE( 12) sEOL( 600) 

TYPE DOUBLE VV eUUsVeX2oXleY2eVY19P1 P2201 sD02sDEL 

NM=0 

JJ=1 © 

READ 1500 eNSsNBRo IOP sHMsAMYK 9 SD oeRI sOMEGAG sOMEGAI » GCDABC eBRC » GRAPH» 
1REM »PMM eHL 

FORMAT(31224F3e0s2F8e0 94F2 e092F4e0) 

IF (NS) 6006002601 

NMBR=NBR+2 

IF (1TOP-1)150121501%1506 

RFAD 15025(SV(1I)sT=1>sNMBR) 

FORMAT ( 22F 30) 

DO 1505 J=1sNMBR 

IF (SV(I9)191503 9150351505 

READ 1504,(SVIP( I oJ) sl #196) 

FORMAT (6F 2-0) 

CONTINUE 

IF (TOP) 1506915061523 

READ 1507,.(SVO(I)»sI=1>sNMBR) 

FORMAT (22F3-0) 

DO 1510 J=l>+NMBR 

IF (SVO(J))15085150821510 

READ 1509,(SVOP(I oJ) sI=196) 

FORMAT (6F 2-0) 

CONTINUE 

IF(NBR) 15139151391511 

READ 1512+5(PHI(1)sIT=1sNBR) 

FORMAT (10F7-3) 

OREAD 15145(0(1)sDI (CE) sRHOC TDs TLE) sTHETACT) os CLX( 1) sCLY( 1) eCLZ(1) oC 
LTX(CT)sCTY( TL) sCTZ( 1) esNNBR( I) sT=19NS) 

FORMAT ( 3F5e292F6e2 s6Fle2sl2) 

READ 1515 » AAMUsAE,AGs00 

FORMAT (F 7e322E10¢e2sF2e0) 

IF (00)15165151621518 

DO 1517 I=1 NS 

AMU (1) =AAMU 

E( 1)=AE 

G(1)=AG 

GO TO 1520 

AMU(1)=AAMU 

E(1)=AE 

G(1)2AG 

READ 15199(AMU(1)sE( 1) sG(I)eT=29NS) 

FORMAT (F7e3s2E10e2) 

IF (OMEGAT) 11219112191521 


62 


000000 
000010 
000020 


~ 000030 


000040 
000050 
000060 
000070 
000080 
000090 
000100 
000110 
000120 
000130 


000140 © 


000150 
000160 
000170 


000180. 


000190 
000200 
000210 
000220 
000230 
000240 
000250 
000260 
000270 
000280 
000290 
000300 
000310 
000320 
000330 
000340 
000350 


~ 000360 


000370 
000380 
000390 


000400 


000410 
000420 
000430 
000440 
000450 
000460 
000470 
000480 
000490 
000500 
000510 
000520 
000530 
000540 
000550 


NAAN 


Yr2t 
1S21 


6091 
6092 


6093 
6094 


6095 


6096 
6097 


6098 


6099 
7000 


7001 
7071 


7072 
7073 


1522 


OMEGAI=.01 
OMEGA(1)=OMFGAI 
AMYK=AMYK=1.40 

OME GAO=OMEGAT 

FFAC=.5 

L=1 

M=] 

MN=0 

AA=0.0 

BB=0.0 

cc=0.0 

FF20.0 

GG=0.0 

RR=0.0 

KK =0 

PP=0, 

KA=0 

KKK=0 * 

AAMU=0.0 

EE20e0 

AL=020 

AD=0.0 

NNN=NBR+1 

IF (PMM#O,1—1.) 609126092 26093 
PmM=-l. 

GO TO 7000 

PM=0. 

GO TO 7000 

IF (PMM#0.01-1e) 609496095 26096 
PM=—1le 

GO TO 7000 

PMzl. 

GO TO 7000 

IF (PMM#0.01—-1e1) 609726098 6099 
PmM=-l. js 

GO TO 7000 

PM=0. 

GO TO 7000 

PMz—-l. 

PPM=PM 

IF(PM) 70739707397001 
IF(REM=1le) 7071970727071 
SS=-l. 

GO TO 7073 

SS=0.6 

DO 1522 [#1lsNS 
OT(I)=01(7T) 

BAMU(T )2AMU(T) 
THETR(IT)STHETACT) 

DIC TI=D(1)-2e0* OI (1) 


INVARIANTS 


DO 108 I=1l»sNS 


IFCCLXCT)+CLYCI)4CLZ 01 +CTXOT) +CTYC TI +CFZ01)) 10851019108 


+ , 


33 


000 560 
000570 
000580 
000590 


* 000600 


000610 


* 000620 


000630 
000640 
000650 — 
000660 
000670 
000680 
000690 
000700 
000710 
000720 


000730 -— 


000740 
000750 
000760 
000770 
000780 
000790 
000800 
000810 
000820 
000830 
000840 


000850 


000860 
000870 
000880 
000890 
000900 
000910 
000920 
000930 
000940 


* 000950 


000960 
000970 
000980 


* 000990 


001000 
001010 
001020 
001030 
001040 
001050 
001060 
001070 
001080 
001090 
001100 
001110 


NAAN 


10] 
102 


103 


104 


10s 


106 
107 


108 


109 


999 


1011 


1021 


AMU( 1) =AMU( 1) #36 1415926548(D(1) ##2- DI(1)**2)/(6912e*32e17#1260) 
IF (RHO(1))102%1039102 

THETA(I1-)=THETA(1) #36 14159265471806. 

TLC I)STHETACI) #ABSF(RHO(I)) 
AJT(1)=301415926548(D(1)#*#4-DI(1)##4)/326 

AJ(I)=AUT(1)/26 

R(LI=LeO/(AJ(I) #EC1)) 

AIX(1)=SORTF((D(1) ##24+DI(1)*#2)/Be) 
ALY(1)=AIX(1)/12414214 

IF (AMYK)10421055104 

Z(1)=120 

IF (RHO(1))105%1069105 

HHM=HM 
HMH=0. 
CALL SUBSEC 
HM=HHM 
IF (NNBR(1))10851075108 

AL=AL+TL(1) 

AD=AD+D(1)#TL(I) 

FE=FE+E(1)*TLO(I) 

AID=AID4+DI(1)*TL(I1) 
AAMU=AAMU+AMU( IT) #TLOI) 

CONTINUE 

AD=AD/AL 

AID=AID/AL 
AFR=36141592654/640*(AD##4-AIDE #4) 
FR=4.73004019##2*#SQORTF (EE ®AFR/ (AAMU#AL £#4) ) 
ADOM=FR/4,.0 

DOM=ADOM 

CR=ADOM#0.001 

BDOM=ADOM 

BR=CR 

ERM=REM 

DO 109 I1=1,5NBR 

PHI (11 )=PHI (11) #36141592654/18000 
IF(1OP-1) 11099995998 

MN=1 


(RHO( 1)» THETACT )oD( 1) oTLO1) sHM9Z(1) sOMEGAT »HMH) 


CONTROL IN PLANE 
P1=1.0 

P2=1.0 

Dl=0.6 

D2=0. 

Q121-0 

Q2=1.0 

Q3=1.0 

B1l=0. 

B2=0-6 

B3=0~6 

FINK=0,. 

IF (SV(1))101191021,.1011 
CALL STAVEC (SV(1)0lsSVI) 
GO TO 1041 

DO 1031 J=1,6 


a4 


001120 
001130 
001140 
001150 
001160 
001170 
001180 
001190 
001200 
001210 
001220 
001230 
001240 
001250 
001260 
001270 
001280 
001290 
001300 
001310 
001320 
001330 
001340 
001350 
001360 
001370 
001380 
001390 
001400 
001410 
001420 
001430 
001440 
001450 
001460 
001470 
001480 
001490 
001500 
001510 
001520 
001530 
001540 
001550 
001560 
001570 
001580 
001590 
001600 
001610 
001620 
001630 
001640 
001650 
001660 
001670 


1031 
1041 
1051 


1061 
1071 
1081 


1085 


1086 
1087 
1000 


WP 


38 
39 
40 
41 
42 
43 
44 


45 


460CALL DISTM 
LISEC(T) saFFACsU) 


SVI(J91l)=SVIP( Jel) 

IF (SVOCNMBR) )10512106191051 © 
CALL STAVEC(SV(NMBR) 01 » SVE) 
GO TO 1081 
DO 1071 J=ls6 

SVE (J91)=SVIP (JU oNMBR ) 

JK=0 

IF (NBR) 1000210001085 

DO 1087 N=2eNNA- 

IF (SV(N))10872108721086 
CALL STAVEC(SV(N) eNo SVIP) 
CONTINUE 

N=2 

T1121 

HH=0-e 

DO 1001 I=1.NS 

MAT=0 

BC=0. 
PCTISCLXCII4CLYCIIV4+CTZ01) 
IF (P(7T) 137238937 

CALL HANGER(CLXCI) eCLYC1) sCTZ(1) VU) 
GO TO 62 

IF (CAMYK) 39941 939 

IF (RHO( 1) 340242940 

IF (Z2(61))43945 943 

IF (RHO( 1) 140248940 

IF (Z(01) 146947246 


CALL CFIELD Ey Sl aE) OMA Y ene ern oNIEN Estee! }'SERee } NOnnenemolesmn Teor 
CALL POINT (AMUCT) sATY C1) sOMEGA(L) o TL 1) 0201) 2R1 9B) 

CALL MATMUL (ABU) 

MAT=1 

GO TO 62 

CALL STIFCO (THETACT ) sRHOCT) 9U) 

GO TO 62 


(AMUCT) sATY (CI) sOMEGA(L) o TL U1) oRI 1) e001) eDIC 1) eGC1)e9SD0eR 


GO TO 62 
CALL RIGID 
GO TO 62 
IF (Z(1)149947 949 
CALL SFIELD (DICT) oTLO IT) oRI1) 0201) 9GC1) »9SDeD( 1) SECT) s9FFACSA) 
GO TO 4&4 

IF (NNBR(IT)) 81080981 

IF (I~-1)63979963 

CALL STATEM(SVI 9UU D1 9D29P19P29PPsPMsSS) 

GO TO 63 

CALL FINMAT(UsUUsV oMAT sA0Z(T)) 

DO 64 J=196 

DO 64 K219e3 

UU (JS oK) =V0UeK) 

IF (PM) 880111384 

IF (FINK) 88284588 

IF (BC) 888801111 

IF (BRC) 88288289 

CALL. BRCOR(R3 e9UU oXXo IK 9 VV 9PM) 

MMM=N-1 


(AMUCT)oTL (1) sOMEGA(L) sATY( 1) oRISU) 


. a” 


BO 


001680 
001690 
001700 
001710 


~ 001720 


001730 
001740 


* 001750 


001760 
001770 
001780 


001790 — 


001800 
001810 
001820 
001830 
001840 
001850 
001860 | 
001870 
001880 
001890 
001900 
001910 
001920 
001930 
001940 
001950 
001960 
001970 


- 001980 


001990 


_002000 


002010 
002020 
002030 
002040 
002050 
002060 
002070 
002080 
002090 
002100 
002110 
002120 
002130 
002140 
002150 
002160 
002170 
002180 
002190 
002200 
002210 
002220 
002230 


NANA 


67 
68 


1001 


998 


1013 


1023 
1033 
1043 
1053 


1063 
1073 
1083 


1093 
1094 


1095 
1002 


DO 90 J=l»s3 

DO 90 K=193 

YY (J9K9MMM) =XX(JoK) 

GO TO 1001 
IF(NNBR(I)-1) 83983265 


CALL STATEB(SVIPsVVeNeYY ol s8RC D1 9D2ePP oPM9SS eG] 9Q2 9Q3 FINK) 
CALL FINBRA(UsVVeVeMAT sAeZ(1I)) 


DO 66 J=196 

DO 66 K=l93 

VV (UeK) RV J9K) 

IF (NNBR(1I))68267967 

IF (NNBR(1)-3)1001968568 
CALL BRANCH(R3sVVePHI(C II) VU) 
N=N+1 

TI2I1I+1 

MAT =O 

BC=le 

GO TO 63 

CONTINUE 

GO TO 2000 


a 


CONTROL OUT OF PLANE 


P1l=1.0 

P2=1.-0 

D1=0O64 

D2=06 

Q1=1.0 

Q2=1.0 

Q3=1.0 

B1l=0.6 

B2=06 

B3=0.6 

FINK=0. 

IF (¢$V0(1))10135102351013 
CALL STAVEO(SVO(1)91sSVI) 
GO TO 1043 

DO 1033 J=1s6 
SVI(J91)=SVOP (Jol) 

IF (SVO(NMBR) )10535106321053 
CALL STAVEO(SVO(NMBR) 91sSVE) 
GO TO 1083 

DO 1073 J=1l+6 

SVE (J91)=SVOP (J »NMBR) 
JK=1 

IF (N8R) 1002100231093 

DO 1095 N=2sNNN 

IF (SVO(N))10959109591094 
CALL STAVEOCSVO(N) »N»oSVOP) 
CONTINUE 

N=2 

Il=1 

HH=0-64 

DO 1003 I=1l»sNS 

BC=06 


836 


002240 
002250 
002260 
002270 
002280 
002290 
002300 
002310 
002320 
002330 
002340 
002350 
002360 
002370 
002380 
002390 
002400 
002410 
002420 
002430 
002440 
002450 
002460 
002470 
002480 
002490 
002500 
002510 
002520 
002530 
002540 
002550 
002560 
002570 
002580 
002590 
002600 
002610 
002620 
002630 
002640 
002650 
002660 
002670 
002680 
002690 
002700 
002710 
002720 
002730 
002740 
002750 
002760 
002770 
002780 
002790 


MAT=0 
P(IT)=CTX(C IV 4CLZ°01)4CTYC1) 
IF(P(IT)) 19201 
1 CALL HANGEOCCTX(1) sTLZ¢(1) sCTY¢(1)9U) 
GO TO 14 
IF (AMYK)391193 
IF (RHO( 1) 45894 
IF(2(01))69596 
CALL STIFOO (THETACT) sRHO(T) 9U) 
GO TO 14 
6OCALL CFIELO(AIJT(I) oRO 199261) 9GC1) 9SOsRHOC(IT) sTHETA(I) sD¢ 1) os DI (1) oF F 
LAC sA) 
7 CALL POINO (AMU(T) sATX(T) SATY (1) sOMEGA(L) oTLOT) 9Z(1) #RI 9B) 
CALL MATMUL (AsBoV) 
MAT=1 
GO TO 14 
8 IF(Z(1))1099910 . 
9 CALL RIGIO (AMUCT) os TLC I) sAIXUI) sOMEGA(L) sAIY(1) sAJT(I) »9RI 9 VU) 
GO TO 14 
LOOCALL DISTMO(CAMU( 1) sAIX(1T) sAITY( I) sOMEGA(L) os TLC I) eAJTII) oR 1) 2D001)90 
LI(T)sGC1)sSOeRI sFFAC 9U) 
GO TO 14 
11 IF CRHOC IT) 1401294 
12 IF (2(1))13999013 
13 CALL SFIELO (DICT) oT LCT) sAJTCI) oROI) 9201) 9GC1) 2SD09D0(1) sFFAC SA) 
GO TO 7 
14 IFCNNBR(I)) 20192009201 
200 IF(I-1) 155202915 
202 CALL STATEM(SVI »UU2D12D29P1sP2,PPsPM,SS) 
GO TO 15 . 
15 CALL FINMAT(UsUUsVeMAT sAsZ(1)) 
DO 16 J=196 
DO 16 K=1>53 
16 UU(JeK) 2V( 59K) 
IF (PM) 20791114985 
1114 IF(FINK) 2079855207 
85 IF(BC) 207920791112 
1112 IF(SRC) 20792072208 
208 CALL BRCOR(R3 sUUsXXe JK 9 VV 9PM) 
MMM=N-1] 
DO 209 J=ls3 
DO 209 K=1,.3 
209 YY(J9K o9MMM) =XX(JeK) 
207 GO TO 1003 
201 IF(NNBR(I)-1) 2049204917 
204 CALL STATEB(SVOP oeVVoNe YY ol BRC» D1sD2sPP, PM»SSQ1» Q22Q3sF INK) 
17 CALL FINBRA(UsVVeVoeMAT sAsZ(1)) 
DO 18 J=196 
DO 18 K=1o93. 
18 VV(JeK) HV IeK) 
IF (NNBR(1))82919919 
19 IFC(NNBR(1)-3) 100358282 
82 CALL BRANCO(R39VVePHI( II) »U) 
N=N+1] 
IT=II+1 
MAT =0 


wm & WA 


oT 


002800 
002810 
002820 


- 002830 


002840 
002850 
002860 
002870 
002880 
002890 
002900 
002910 
002920 
002930 
002940 


002950 | 


002960 
002970 
002980 - 
002990 
003000 
003010 
003020 
003030 
003040 
003050 
003060 
003070 
003080 
003090 
003100 
003110 


-003120 


003130 
003140 
003150 
003160 
003170 
003180 
003190 
003200 
003210 
003220 
003230 
003240 
003250 
003260 
003270 
003280 
003290 
003300 
003310 
003320 
003330 
003340 
003350 


C 
c 
.S, 


1007 


2000 CALL DELMA(UUsSVEeVeX2 sV2eKKKoX1 oF F »PP »sHHeRReSSeVl sREMoFINK oPMo D1 » 
102) 


302 


SSODEL(LI=AVE Le lI #V6292)*V(303)-VE1L eS FVI2 022 HV (391 +V0192)*#V02 03) #V03 
LelI-“V(192)*V06 291) *V0303)4+V06 103) V0 302) FV0 20) P-VE lel) FV03902) FV(2 93) 


21 


2 900MEGAM(M) =OMEGAC+ ( OMEGA (L) -OMEGAC) *ABSF (DELC)/ (ABSF (DELC) + 
LABSF(DEL(L))) 


30 


31 


73 


86 
32 


BC=le 
GO TO 15 
CONT INUF 


ITERATION 


IF (FINK) 87287586 

IF (PM) 99299998 

IF (KKK) 7030270302301 
GO TO 6001 ; 

LF (FF-12e¢0) 7051270607051 
KK=] 

M=M—-] 

GO TO 3000 

IF (HH-12¢0) 300260015300 
IF (SS) 30343022302 

DEL (L)=X2-Y2, 

D1=X1 

D2=(X24+Y2)/2e0 

GO TO 61 

DEL(L)=X1-Y1 
01=(X14+Y1)/2e0 

D2=X2 

GO TO 6} 


GO TO 61 
PDEL=DEL(L) 

GO TO 25 

IF (AA) 23923924 
PDOEL=DEL(L-1) 

GO TO 25 
PDEL=DELC 

IF (PDEL#DEL(L))26226927 
IF (AA) 28928957 
BB=0 

IF (AA) 359355952 
AA=1,.0 

BB=1.0 
DELC=DEL(L-1) 
OMEGAC =OMEGA(L=-1) 
GO TO §9 


GO TO 32 


OMEGA(L+1)=OMEGA(L)-029*#DOM 


DOM=0.5*DOM 

L=L+1 

BL=L 

IF(HL-BL) 705397053586 
KA=1 ; 

GO TO 3000 

IF (JK) 1000+ 100021002 
BM=M 


003360 
003370 
003380 
003390 
003400 
003410 
003420 
003430 
003440 
003450 
003460 
003470 
003480 
003490 
003500 
003510 
003520 
003530 
003540 
003550 
003560 
003570 
003580 
003590 
003600 
003610 
003620 
003630 
003640 
003650 
003660 
003670 
003680 
003690 
003700 
003710 
003720 
003730 
003740 
003750 
003760 
003770 
003780 
003790 
003800 
003810 
003820 
003830 
003840 
003850 
003860 
003870 
003880 
003890 
003900 
003910 


 —~—* = 


33 


34 
340 


35 
50 


a1 
52 
53 


54 


5 SOOMEGAM( M) =OMEGA(L )+(OMEGA(L- Ly =OMEGACL 1) SABSF UREA aly deliiiabali alee 
1)+ABSF(DEL(L=-1))) 


3000 


401 


IF (HM=BM) 3000 23000 933 
M=M+] 
IF (M-3) 34090349 34 


ADOM= ( OMEGAM( M~1)-~OMEGAM( M- ~21)*0. ] 


AA=020 
BB=0 
GC=0 
GG=0.0 
DOM=ADOM 
DELC=0 


OMEGA (L+1 ) =OMEGAM( M- 1)+2e*CR 


OME GAO=OMEGA(L+1) 
GO TO 31 


IF (OMEGA(L ) -OMEGAG) 53260260 . 
OMEGA(L+1)=OMEGA(L )-(OMEGA(L )-OMEGAC)/220 


IF COMEGA(L ) -OMEGA(L+1)-CR)51 951 031 


CC=160 

GO TO 31 

DELC=DEL(L) 
OMEGAC=OMEGA(L) 
OMEGA(L+1)=OMEGA(L)+DOM 
GO TO 32 

IF (DELC*#DEL(L) 129929955 


GO TO 32 

IF (OMEGA(L ) -OMEGAO) 21921922 

IF (8B)58258,50 

BB=1.0 

IF (OMEGA(L )-OMEGA(L=-1) -CR) 29929930 
M=M-] 

GO TO 3000 

IF (6G) 70970072 

IF(DEL(LII 77971977 

GG=1.0 


OMEGA(L+1)=OMEGA(L)-CR 
GO TO 31 

IF (GG-120) 73973974 
GG=2e0 

HDEL=DEL(L) 
OMEGA(L+1)=OMEGA(L)+2e¢#CR 
GO TO 31 
TOEL=DEL(L)*HDEL 

IF (TDEL) 75976076 
OMEGAM(M) =OMEGA(L=-2) 
GO TO 33 

KK=1 

M=M-1] 

GO TO 3000 

IF (CC)56256954 


OUTPUT 
TF (NM) 400 0400 2403 


PRINT 401 
FORMAT (1H1//) 


87 


003920 
003930 
003940 
003950 
003960 
003970 
003980 
003990 
004000 
004010 
004020 
004030 
004040 
004050 
004060 
004070 
004080 
004090 
004100 
004110 
004120 
004130 
004140 
004150 
004160 
004170 
004180 
004190 
004200 
004210 
004220 
004230 
004240 
004250 


004260 


004270 
004280 
004290 
004300 
004310 
004320 
004330 
004340 
004350 
004360 
004370 
004380 
004390 
004400 
004410 
004420 
004430 
004440 
004450 
004460 
004470 


PRINT 402 
402 FORMAT (15H PROGRAM VIPIPEs35X» 15HYe Se KIM NHA-3938X» 13H MARCH 


1 1966//910X»2108H MODAL FREQUENCIES OF A PLANAR PIPING 
2 SYSTEM ARE DETERMINED BY AN ITERATIVE PROCEDURE USING THE /29H ME 
3THOD OF TRANSFER MATRICES e///119H SH HEHEHE EEEE RHE HE EERE EE HEE EEE 


y Hea AE AE SE A aE EE da aE aE aa a ae aa aa aa aa 


> hela Miet ated R elie teMetetMetMelielial 
IF (NM) 4058» 4058 2404 
PRINT 405 
FORMAT (1H1//) 
NM= 1] 
4059 IF( JK) 4060+ 406024065 
4060 PRINT 40709JJ 
4OTOOFORMAT (50X»8HPROBLEM 129//35Xs47HIN PLANE MODE FREQUENCIES 
INS PER SECOND) //45X e4HMODE 5 15Xe10H FREQUENCY /) 
GO TO 4079 
4065 PRINT 40759JJ 
40750FORMAT(50X»8HPROBLEM I120//33Xe51HOUT OF PLANE MODE FREQUENCIES 
1DIANS PER SECOND) //45Xs4HMODE 915X910H FREQUENCY /) 
MH=M 
PRINT 4089(M»sOMEGAM(M) »M=1] »MH) 
FORMAT (46X912917X9F14e8 /) 
IF(KA) 5900,5090,5010 
PRINT 5020 
FORMAT(10Xe53HPROBLEM TERMINATFD DUE TO OMEGA STORAGE LIMITATION. 
17 ) 
5000 IF(KK) 4080» 409054080 
4080 PRINT 4085 
4O850FORMAT(10X»69HPROBLEM TERMINATED DUE TO SIGNIFICANT FIGURE LIMITAT 
LION OF COMPUTER. /) 
4090 IFCAMYK)412 54095412 
409 PRINT 411 
4110FORMAT (16X»53H 
1ND ) 
GO TO 414 
PRINT 413 
FORMAT (24X%»s45HA DISTRIBUTED MASS APPROACH WAS EMPLOYED AND ) 
1F(S0)4155415 9420 
415 IF(RI)41694169418 
416 PRINT 417 
417OFORMAT (2Xe67HTHE EFFECTS OF SHEAR DEFLECTION AND ROTARY INERTIA W 
1ERE CONSIDERED. /) 
GO TO 4610 
418 PRINT 419 . 
4190FORMAT (3Xe63HTHE EFFECTS OF SHEAR DEFLECTION WERE CONSIDERED WHIL 
1E THOSE OF /3Xs35HROTATIONAL INERTIA WERE NEGLECTED. /) 
GO TO 4610 
420 IF(RI)42194212423 
421 PRINT 422 


//) 
403 
404 
405 

4058 


(RADIA 


(RA 
4079 
408 


5010 
5020 


A LUMPED MASS APPROACH WAS EMPLOYED» A 


412 
413 
414 


42Z220FORMAT (3Xe972HTHE EFFECTS OF ROTATIONAL INERTIA WERE CONSIDERED WH 
LILE THOSE OF SHEAR /3Xs27HDEFLECTION WERE NEGLECTED. /) 
GO TO 4610 


423 PRINT 424 
4240FORMAT( 3X966HTHE EFFECTS OF SHEAR DEFLECTION AND ROTARY INERTIA WE 
1RE NEGLECTED. /) 
4610 IF(PM) 425942594713 


ate f 


70 


004480 
004490 
004500 
004510 


. 004520 


004530 
004540 
004550 
004560 
004570 
004580 
004590 
004600 
004610 
004620 
004630 
004640 
004650 
004660 
004670 
004680 


- 004690 


004700 
004710 
004720 
004730 
004740 
004750 
004760 
004770 
004780 
004790 
004800 
004810 
004820 
004830 
004840 
004850 
004860 
004870 
004880 
004890 
004900 
004910 
004920 
004930 
004940 
004950 
004960 
004970 
004980 
004990 
005000 
005010 
005020 
005030 


4713 
4611 
4613 
4623 


IF (SS) 4611946129%612 

IF (PP=1le) 461329461494615 

PRINT 4623 

FORMAT(24Xe935HPESTEL MAHRENHOL TZ 
GO TO 425 

PRINT 4624 

FORMAT (24Xs35HPESTEL MAHRENHOL TZ 
GO TO 425 

PRINT 4625 
FORMAT (24Xe35HPESTEL MAHRENHOL TZ 
GO TO 425 
it SO 1 6s) 
PRINT 4626 
FORMAT ( 24Xs35HPESTEL MAHRENHOL TZ 
GO TO 425 

PRINT 4627 

FORMAT( 24Xs35HPESTEL MAHRENHOL TZ 
GO TO 425 

PRINT 4628 

FORMAT (24X935HPESTEL MAHRENHOL TZ 
IF (PM) 90102*9020*9000 

PRINT 9012 

FORMAT (24Xs20HDELTA METHOD USEDe /) 

GO TO 9000 

PRINT 9022 

FORMAT(24Xs29HMODIFIED DELTA MFTHOD USED. /) 
PRINT 440 


METHOD(A) USEDe /») 
4614 


4624 METHOD(B) USED. /) 
4615 


4625 METHOD(C) USEDe /) 
4612 
4616 


4626 


4616946174618 


METHOD(D) USEDe pie 
4617 | 


4627 METHOD(E) USEDe /) 
4618 
4628 

425 
9010 


9012 


METHOD(F) USEDe /). 


9020 
9022 
9000 


4400FORMAT ( 118H HEE EEE EEE REE EEEE EE EERE EERE EEE EEE EEE EEE EEE. 


1 SEEEEEEEEEEEEEEEE EEE EEE EEE ERE REESE REESE EERE ERE EE REE EH REE EE ) 


IF (GDABC=-1¢. )45112452094401 
4401 PRINT 4402 


44020FORMAT(35H SECTION LENGTH AND SUBSECTION DATA/23X914HSECTION NUMBE 


1ReBXe25HLENGTH OF SECTION( INCHES) »5Xe21HNUMBER OF SECTIONS = 
PRINT 443 ¢( ToTL( IT) oZ(1)eT21leNS ) oe 
443 FORMAT (29Xe13e23XeFle2921X9F4,0) 
LM=L 
IF (PM) 7054,7054,.7055 
7055 PRINT 44] : 
4410FORMAT (11H GRAPH DATA /22Xs16HITERATION NUMBER »4Xs30HFREQUENCY (R 
LADIANS PER SECOND) 93X230H VALUE OF REMAINDER ) 
GO TO 7057 
7054 PRINT 7056 
TOS60FORMAT (11H GRAPH DATA /22X+16HITERATION NUMBER »s4Xs30HFREQUENCY (R 
LADIANS PER SECOND) 93Xs30HVALUE OF FREQUENCY DETERMINANT ) 
PRINT 4449(L»sOMEGA(L) »sDEL(L) sLel»LM) 
FORMAT (29X913923Xe9F1l1le6918X9E15¢e8) | 
PRINT 442 
FORMAT (1Xs10HINPUT 
PRINT 4521 
45210FORMAT (1X s9HCOMPONENT 9 4X s8HDIAMETERs5Xs14HWALL THICKNESS »5X s6HLEN 
1GTHsS5Xes9HRADIUS OF s9Xs8HINCLUDED 9 5X»s THDENSITY 9 5Xs THELAST IC 95Xe5HSH 
Z2EAR /S7TXe9HCURVATURE oe 5Xe12HANGLE OF ARC+16Xs 7THMODULUS »5X » THMODULUS 
3 /) 
DO 4525 I=1l>sNS 
IF (P (1) )45249452454523 
4524 PRINT 452291 eD( 1) eDT( I) eTL( I) sRHOC I) sTHETR(I) sBAMU(I)SE(I) sG(1) 


7057 
444 
4520 


442 DATA /) 


Z| 


005040 
005050 
005060 
005070 


005080 


005090 
005100 
005110 
005120 
005130 
005140 
005150 
005160 
005170 
005180 
005190 
005200 
005210 
005220 
005230 
005240 
005250 
005260 
005270 
005280 
005290 
005300 
005310 
005320 
005330 
005340 
005350 
005360 


_ 005370 


005380 
005390 
005400 
005410 
005420 
005430 
005440 
005450 
005460 
005470 
005480 
005490 
005500 
005510 
005520 
005530 
005540 
005550 
005560 
005570 
005580 
005590 


45220FORMAT (3X91398XoF 76 45 NOGeoe 39 TX9F Bed o4Xs F8e399XeFbe2e8XoF7e294Xs 005600 


Fa 


1E80294X9E8e2) 005610 

GO TO 4525 005620 

4523 SH=le : 005630 
4525 CONTINUE 005640 
IF (SH) 453094530 94526 . 005650 

4526 PRINT 4532 005660 
4532 FORMAT (//1X»7HHANGERS /) 005670 
PRINT 4527 005680 

452 70FORMAT (2Xs9HCOMPONENT 914X 9 3HCLX914X 9 3HCLY 914X93HCLZ914X93HCTX914X 005690 
1s3HCTY914X93HCTZ /) 005700 
4528 DO 4530 I=1l9NS -_ 005710 
IF (P(1) 1453094530 94529 005720 

4529 PRINT 453loTeCLX( 1) sCLY(I) sCLZ(1) oCTX( 1) sCTY(1) sCTZ¢1) 005730 
4531 FORMAT(3X913d15X oF 1lO0e2e7TXoFlOe2 os TXeFlOe2eTXeFlOe2oTXeFlOe2eTX9F1006e 005740 
ze) 005750 
4530 CONTINUE 005760 
PRINT 448 005770 

448 FORMAT (//20H BOUNDARY CONDITIONS ) 005780 
PRINT 449e9(SVI(191)91=196) 005790 

449 FORMAT (5H SVI=6F3e0) 005800 
PRINT 45le(SVE(1901)s7T=196) 005810 

451 FORMAT (5H SVE=6F3.0) 005820 
IF (NBR) 4511945114502 005830 

4502 IF (JK)450554505 94515 005840 
4505 PRINT 4510e((SVIP(19J) 91=196) »sJ=29NNN) 005850 
4510 FORMAT (6F3-0) 005860 
GO TO 4511 005870 

4515 PRINT 45179 ((SVOP(1sJ)91=196) »J=29NNN) 005880 
4517 FORMAT (6F3.00) 005890 
4511 CONTINUE 005900 
SORT -005910 

IF (GRAPH=1e) 6001600026001 005920 

6000 LM1=LM=~1 005930 
DO 6100 L=1+sLM 005940 
EDL(L)=DEL(L) 005950 

6100 CONTINUE 005960 
DO 6600 L=1,LM) 005970 
LPl=L+] 005980 

DO 6600 M=LP] LM 005990 

IF (OMEGA(L) -OMEGA(M)) 660056600 »6003 006900 

6003 TEMP=OMEGA(L) | 006010 
OMEGA(L)=OMEGA(M) 006020 
OMEGA(M)=TEMP 006030 
TEMP=EDL(L) 006040 
EDL(L)=EDL(M) 006050 
EDL(M)=TEMP 006060 

6600 CONTINUE 006070 
C GRAPH 006080 
NUMPTS=60 006090 

K=0 006100 

DO 6006 J=1912 006110 

6006 ITITLE(J)=8H 006120 
ITITLE(3)=8H KIM Y S 006130 
ITITLE(8)=8HFREQUENC 006140 
ITITLE(9)=8HY VS R 006150 


ITITLE(10)=8HEMAINDER : “ 006160 


LT1=LM/NUMPTS 006170 
IFC(LTL) 602426025 26024 : 006180 

6025 NUMPTS=LM 006190 
LABEL=4H 006200 

L=1 . 006210 

CALL DRAW (NUMPTSeOMEGA sEDL 10 s0eLABEL so ITITLE 0909090909009» 006220 
110sKeLAST) 006230 

GO TO 6001 . : * 006240 

6024 DO 6020 J=lel? -: , 006250 
LT2=0 006260 
LTZ2=LT2+J*NUMPTS 006270 

IF (LT2-LM) 602026021 96022 006280 

6020 CONTINUE . . . ii 006290 
6021 LT=LT1 i: . 006300 
GO TO 6023 ° 006310 

6022 LT=LT1+1 006320 
NUMEND=LM-LT1*NUMPTS 006330 

6023 L=l . 006340 
DO 7021 JT=lsLT 006350 . 

IF (I-1) 6008 2604126008 , 006360 

6008 IF(I-LT) 8011560116011 | 006370 
BO11 IF(1-3) 6042604398012 | 006380 
8012 IF{ 1-5) 6044604558013 | " 006390. 
9013 IF(I-7) 6046 2604798014 ; 006400 
8014 IF(I-9) 6048 9604928015 006410 
8015 IFCI-LT) 6050260116001 006420 
6041 LABEL=4H ONE 006430 
GO TO 6009 006440 

6042 LABEL=4H TWO 006450 
GO TO 6010 006460 

6043 LABEL=4HTHRE 006470 
GO TO 6010 006480 

6044 LABEL=4HFOUR . 006490 
GO TO 6010 006500 

6045 LABELZ4HFIVE 006510 
GO TO 6010 006520 

6046 LABEL=4H SIX 006530 
GO TO 6010 006540 

6047 LABEL=4HSEVN 006550 
GO TO 6010 006560 

6048 LABEL=4H EIT 006570 
GO TO 6010 006580 

6049 LABEL=4HNINE 006590 
GO TO 6010 006600 

6050 LABEL=4H TEN 006610 
GO TO 6010 006620 

6009 DO 6031 N=zleNUMPTS 006630 
FDL(N)=EDL(L4+N—-1) 006640 

6031 OMEGA(N) =OMEGA(L4+N-1) 006650 
PRINT 6027»sNUMPTSsEDL(1) 006660 

6027 FORMAT(1H10///////5X THNUMPTS= 9 1395Xe7THEDL(1)= 961209) 006670 
CALL DRAW (NUMPTS » OMEGA sEDL 9020,LABEL sITITLE%0 909090009009 006680 
110sKsLAST) 006690 

GO TO 7021 006700 

6010 LeL+NUMPTS 006710 

pat” 


We 


OMEGAX=OMEGA(L) 
DO 6030 N=leNUMPTS 
EDL (N)=EDL(N+L—-1) 
OMEGA(N) =OMEGA(L+N-1) 
6030 OMEGA(N)=OMEGA(N) -OMEGAX 
PRINT 6032+NUMPTS»EDL(1) 
6032 FORMAT(///////5X97THNUMPTS=+ 13 »95Xs7HEDL(1)=eE1209) 
CALL DRAW(NUMPTS OMEGA sEDL 00 290 sLABELsITITLE 2020909209090 99» 
110K »sLAST) 
GO TO 7021 
6011 L=L+NUMPTS 
OMEGAX=OMEGA(L) 
DO 6033 N=1l»NUMEND 
EDL (N) =EDL(N4+L—1) 
OMEGA(N) =OMEGA(N+L—1) 
6033 OMEGA(N) =OMEGA(N) -OMEGAX 
PRINT 6051»sNUMEND,EDL(1) 
6051 FORMAT(///////5XsTHNUMENDE 9139 5Xs7THEDL(1)*2E12e9) 
LABEL=4HLAST 
CALL DRAW(NUMEND»s OMEGA s EDL 2020 sLABELsITITLE 20202020 209099>% 
1109K esLAST) 
GO TO 6001 
7021 CONTINUE 
6001 IF(HH-1-0) 802458023 58023 
8023 MH=M 
LM=L 
8024 DO 446 M=1>+MH 
446 OMEGAM(M) =020 
DO 447 L2leLM 
OMEGA(L)=020 
EDL (L)#0. 
447 DEL(L)20.0 
DELC=020 
OMEGAC20.0 
IF(HH-12e0) 8001280228001 
8001 IF(PM) 8031580318033 
8033 IF(REM) 8031%803258032 
8031 IF(MN) 802258022 9452 
8032 IF(REM=1.) 8041280428043 
8041 IF(PP-1) 8043280438072 
8072 IF(SS) 8073580748074 
8073 SS*0-6 
GO TO 8022 
8074 IF(SS=1e) 804228042 8075 
8075 SS=-le 
PP=#0. 
REM=-1l, 
GO TO 8033 
8042 SS=SS-le 
IF(SS) 8051,58052,8053 
8051 SS=le 
GO TO 8022 
8052 SS=26 
GO TO 8022 sd 
8053 SS=-le 
REM=—le 


94 


006720 
006730 
006740 
006750 
006760 
006770 


. 006780 


006790 
006800 
006810 
006820 
006830 
006840 
006850 
006860 
006870 
006880 
006890 
006900 
006910 
006920 
006930 


— 006940 


006950 
006960 
006970 
006980 
006990 
007000 
007010 
007020 
007030 


‘007040 


007050 
007060 
007070 
007080 
007090 
007100 
007110 
007120 
007130 
007140 
007150 
007160 
007170 
007180 
007190 
007200 
007210 
007220 
007230 
007240 
007250 
007260 
007270 


8043 
8061 
8062 


8063 


452 


8022 


7040 
8191 
8192 
8091 
8081 


8092 
8082 


8093 
8084 


8291 
8281 


8292 
8282 


8293 
8283 
8021 
8035 
8036 
8037 

453 
6081 
6084 


GO TO 8033 

PP=PP-1. 

IF (PP) 8061580628063 
PP=le 

GO TO 8022 

PP=z=2e 

GO TO 8022 

PP=06 

REM=-le 

GO TO 8033 

PP=0. 

RR=0. 

SS=-le 

L=] 

Mz] 

AA=0.0 

BB=0.0 

CC=040 

FF=0.0 

GG=0-0 

KK=0 

KA=0 

KKK=0 

ADOM=BDOM 

CR=BR 

OMEGA(1)=OMEGAITI 
OMEGAO#OMEGAIT 

DOM=ADOM 

IF (HH=-12¢0) 80219704098021 
IF (SS) 8191581928192 
IF(PP=le) 8091809248093 
IF(SS—le) 82919829258293 
PRINT 8081 

FORMAT (1H19////18HMETHOD 
GO TO 8021 

PRINT 8082 

FORMAT (1H19////18HMETHOD 
GO TO 8021 

PRINT 8084 

FORMAT (1H19////18HMETHOD 
GO TO 8021 

PRINT 8281 

FORMAT (1H19////18HMETHOD 
GO TO 8021 

PRINT 8282 

FORMAT (1H19////18HMETHOD 
GO TO 8021 

PRINT 8283 

FORMAT (1H19////18HMETHOD 
IF (PM) 453545358035 
IF(REM) 80379803698036 
IFC JK) 99899999998 
REM=ERM 

IF(PM) 6081560826080 

IF (PMM-1.) 608456080 26084 


FAILED. 


FAILEDe 


FAILEDe 


FAILED. 


FAILEDs 


FAILED. 


IF(PMM#00¢1l-100¢1) 608626085 96086 


95 


4//) 


4/7) 


///y) 


4/1) 


4/7) 


//y 


007280 
007290 
007300 
007310 


- 007320 


007330 
007340 
007350 
007360 
007370 
007380 
007390 
007400 
007410 
007420 
007430 
007440 
007450 
007460 
007470 
007480 
007490 
007500 
007510 
007520 
007530 
007540 
007550 
007560 
007570 
007580 
007590 
007600 


-007610 


007620 
007630 
007640 
007650 
007660 
007670 
007680 
007690 
007700 
007710 
007720 
007730 
007740 
007750 
007760 
007770 
007780 
007790 
007800 
007810 
007820 
007830 


30 


PM=le 
GO TO 8036 
PM=06 
GO TO 8036 


IF (PMM#O01-1el) 6080%6080 +6085 


IF (MN) 9001590018034 
MN=0 

PM=PPM 

GO TO 998 

JJ=eJJ+1 

DO 454 I=1 NS 

P(1)=0O6 

Z(1)=020 

GO TO 450 

END 


X 


SUBROUTINE SUBSEC (RHO sTHETA sD» TL eHMeZ »OMEGA I »HMH) 


RHOD#ABSF (RHO) 

IF (HMH) 41240041 

IF (HM—-520)89807 
IF{OMEGAI-800.) 109797 
HM=6. 

GO TO 10 

HM2HMH 

RA=TL/D 

IF (RHOD-02) 39394 

IF (RA=10¢)59596 

Z=06¢ 

RETURN 

IF(RA=30e) 15915916 
Z=220 

RETURN 

IF (RA—-60)17917918 
22340 

RETURN 

IF ({HM-30) 19192 

2=620 

RETURN 

Z=2e¢0*HM 

RETURN 

RD=RHOD/D 

IF (RO-10)19919923 

Z=0 

RETURN 

IF (RD-30) 24924925 

IF (THETA 044) 26026027 
Z=06 

RETURN 

IF ( THETA~12¢57) 28928929 
Z=420 

RETURN 

2=620 

RETURN 

IF (RD-60) 30930931 

IF (THETA 044) 32932933 


76 


007840 
007850 
007860 
007870 
007880 
007890 
007900 


- 007910 


007920 
007930 
007940 
007950 
007960 
007970 
007980 
007990 
008000 
008010 
008020 
008030 
008040 
008050 
008060 
008070 
008080 
008090 
008100 
008110 
008120 
008130 
008140 
008150 
008160 


‘008170 


008180 
008190 
008200 
008210 
008220 
008230 
008240 
008250 
008260 
008270 
008280 
008290 
008300 
008310 
008320 
008330 
008340 
008 350 
008360 
008370 
008380 
008390 


BZ 


32 
34 


Eye) 


BP 
36 


& WP 


ow 


2=4.29 

RETURN 

IF ( THE TA-1-657)34934, 35 
276.0 

RE TURN 

2=8.0 

RETURN 

IF (THE TA=244) 369 3693 
224.20 

RETURN 

AAD 


SUBROUTINE DISTM (CAMUsATYsOMEGAsTL»ReDeDI sGeSDoRI sEsFFACsU) 


DIMENSION U(606) 
ARA=(32141592654/4 .0)#* (D¥*¥2-DI #*2) 
AREA=ARA#FFAC 
BE=TL*¥OMEGA#SQORTF (AMU/ (E*#ARA) ) 
P&4=R¥EAMUFTL ¥#4 #OMEGA#¥2 
IF(RT)Lols2 
T=R*¥AMU*®(AIY#OMEGA*#TL) ®*2 

GO TO 3 

T=0.-0 

IF(SD)40495 

S= (AMU* (OMEGA*TL) *##*2)/ (G*#AREA) 
GO TO 6 

$=0.0 

SPT=S+T 
VV=SORTF(P4+(S-T)**#2/4e0) 
AL1=SORTF(VV--5*SPT) 

AL2=SORTF (VV+e5*SPT ) 
CL=LleO/( ALL ¥*#¥2+AL2** 2) 
CHLI=(EXPFCAL1)+1e¢/EXPFI(AL1))/26¢ 
SHLL=CEXPFI(ALL)I—Le/EXPFIAL1))/27- 
CL2=COSF(AL2) 

SL2=SINFCAL2) 

SBE=SINF(BE) 

CBE=COSF(BE) 

COO=CL¥*¥ (AL 2##2*CHL1+AL1*#2*CL2 ) 
COL=CL¥ ( (AL2**2*SHL1)/AL1+(AL1#*2*SL2)/AL2) 
CO2=CL¥* (CHL1-CL2) 
CO3=CL*(SHL1/AL1-SL2/AL2) 

NO 7F J=196 

DO 7 K=196 

U(JoK) 2000 

Ut(1l»s1)2CBE 

U(1s6)=TL*SBE/ (BE*#ARA*E ) 
U(292)=CO0-S*CO02 
U(293)=TL*(CO1-SPT*CO3) 
U(294)=R*ECO2ZRTL#E¥2Z 

U(2e5) =(R*¥TL#¥¥3/P4)*( (P44+S**2)*#CO3-S*CO1 ) 
U(392)=P4*CO3/TL 

U( 393) =CON0-T#CO2 
U(394)=TL#R*¥(CO]1-T*¥C0O3) 
U(395)=U(294) 

RR=1-0N/R 


7a 


008400 
008410 


“ 008420 


008430 
008440 
008450 
008460 


008470 


008480 
008490 


* 008500 


008510 
008520 
008530 
008540 
008550 
008560 
008570 
008580 
008590 
008600 
008610 
008620 | 
008630 
008640 
008650 
008660 
008670 
008680 
008690 
008700 
008710 
008720 


. 008730 


008740 
008750 


“008760 


008770 
008780 
008790 
008800 
008810 
008820 
008830 
008840 
008850 
008 860 
008870 
008880 
008890 
008900 
008910 
008920 
008930 
008940 
008950 


ne w 


U(492)=P4#*RR*¥CO2/TL#¥*2 
U(493)=RR¥((P44+T#*2) *CO3-T#CO)])/TL 
U(494) =U( 393) 

U(455)=U(293) ; 
U(592)=P4#*RR*(CO1-S*CO3)/TL¥*3 
U(593)=U(492) 
U(594)=U(392) 

U(595)=Ul( 292) 

U(69T)=-AMU*TL *#OMEGA* #2 *SBE/BE 
U(6s6)=CBE 
RETURN 

END 


SUBROUTINE DISTMO (AMUsAIX sALY sOMEGAs TL eAIUT eReDeDI »GeSDoRI oF FAC oU) 


DIMENSION U(606) 

W=le/(G#AJT). 
AREA=(30141592654/4.0) *(D¥#2-DI] #¥2) #F FAC 
BE=+SQRTF(AMU * ( T L *AIX*OMEGA) ##2#W) 
P&=R*AMU * TL ¥#4*OMEGA*##2 
IF(RI) 29293 
T=R¥AMU*#(AILY*OMEGA*TL ) ##2 

GO TO 4 

T=0.0 

IF (S0)59596 

S=(AMU#* (OMEGA*TL) ##2)/ (G#AREA ) 

GO TO 7 

$20.0 

VV=+SORTF (P44(S-T) ##2/420) 
SPT=S+T 

AL1=+SQRTF(VV-.5*SPT) 
AL2=+SQRTF(VV+e5*SPT) 
CL=1e¢0/(AL1*#2+AL2** 2) 
CHLIF(EXPFI(AL1)+1le/EXPF(AL1))/26¢ 
SHL1I=(EXPFU(AL1)~1le /EXPFI(AL1))/26 
CL2=COSF(AL2) 

SL2=SINF(AL2) ~ 

SBE=SINF (BE) 

CBE=COSF( BE) 

COOFCL#(AL2*#2*CHL1] + AL1**2*CL2) 
CO]lSCL#( (AL2**2*SHL1)/AL1 + (AL1**2*SL2)/AL2) 
CO2=CL*#(CHL1I-CL2) 
CO3=CL*¥(SHLI/AL1-SL2/AL2) 

DO 1 J=1s6 

DO 1 K#1le6 

U(JsK) 2020 

U(1ls1)=CBE 

U(ls2)=TL#*W*SBE/BE 
U(2e1l)=-AMU * TL *#AIX##2*SBF/BE 
U(292)2CBE 

U(353)=CO0-S#*CO02 

U(304)= TL *(CO1-SPT#CO3) 
U(395)=R*¥*CO2*TL##2 
U(36)2(R*TL##*3/P4)*#(-S#CO14+(P44+S5*#*2)*C03) 
U(493)=P4*CO3/TL 

U(4s4)=CO0-T#*CO2 


98 


008960 
008970 


008980 
008990 
009000 
009010 
009020 
009030 
009040 
009050 
009060 


- 009070 


009080 
009090 
009100 
009110 
009120 
009130 
009140 
009150 
009160 
009170 
009180 
009190 
009200 
009210 
009220 
009230 
009240 
009250 
009260 
009270 
009280 
009290 
009300 
009310 
009320 


009330 


009340 
009350 
009360 
009370 
009380 
009390 
009400 
009410 
009420 
009430 
009440 
009450 
009460 
009470 
009480 
009490 
009500 
009510 


U(455)= TL*R*(CO1-T*CO3) 
U(496)=U(395) 

RR=1e2e0/R 
U(593)=P4*RR*¥CO2/TL¥*¥2 

U( 594) =RR#(-THCO1+ (P4FT#H2 )*CO3)/TL 
U(595)=U(494) 

U(5 96) =U(354) 

U( 693) =P4*RR*(COl- ~S*CO3)/TL¥*3 
U(694)=U(503) 

U(695)=U(4_93) 

U(696)=U( 393) 

RETURN 

END 


SUBROUTINE RIGID (AMUeTLe OMEGA ALY oRI 9 U) 
DIMENSION U(696) 

AM=AMU*TL 

DO 1 J=1l6 

DO 1 K=196 

Ul JoK)=020 

U(1lel)=1-20 
U(691)=-AM*OMEGA**2 
U(606)=1e9 

U(292)=120 

U(293)=TL 

U(393)=1.20 

U( 492) <AM*#¥TL#¥OMEGA ** 2/2 60 
IF(RI)29293 

U( 493) =AM*OMEGA#*#2*( TL¥*2/6e0-ALY**2) 
GO TO 4 . 

U( 493) =AM* (OMEGA#TL) **27600 
U(494)=1,0 

U(4,5)=TL . 
U(502)2-U(651) 
U(593)=U(492) 

U(5»95)=1e0 

RETURN 

END 


SUBROUTINE RIGIO (AMUsTLsAIX OMEGA PATY sAJTsSDRI 9 U) 


DIMENSION U(606) 

DO 2 J=ls6 

DO 2 K=196 

Ul JoK)=020 

AM=AMU* TL 

U(1ls1)2#1-20 

U( 291) =-AMU*TL¥(AITX*OMEGA) ¥#2 
U(292)=1¢20 

U(303)=1.20 

U(304)=2TL 

U(404)2=1.0 

U( 593) =AM*TL*¥OMEGA*#* 2/2 20 
IF(SORI-O.) 30230931 


30 U5 04) =AM#OMEGA##2*( TLE¥2/600-ALYE*2 ) 


a 


009520 


009530 
009540 
009450 
009560 
009570 
009580 


. 009590 


009600 
009610 
009620 
009630 
009640 
009650 
009660 
009670 
009680 
009690 
009700 
009710 
009720 
009730 
009740 
009750 
009760 
009770 
009780 
009790 
009800 
009810 
009820 
009830 
009840 


-009850 


009860 
009870 
009880 
009 890 
009900 
009910 
009920 
009930 
009946 
009950 
009960 
009970 
009980 
009990 
010000 
010010 
010020 
010030 
010040 
010050 
010060 
010070 


GO TO l 


U( 594) =AM*OMEGA**2*(TL#¥*2/60) 


10555) =1.0 
U(5e6)=TL 


U( 653) =AM*OMEGA** 2 


U(654)=U(503) 
U(696)=1¢40 
RETURN 

END 


SUBROUTINE STIFCO 
DIMENSION U(696) 


DO 1 J=1l6 
DO 1 Ktl»e6 
U(JIeK)=020 
IF (RHO) 29293 


V=1.0 
CT=COSF(THETA) 
ST=SINF(THETA) 
Utlel)=CT 
Ut(le2)2-V*¥ST 
U(201)2V*ST 

Ui 2e2)=CT 
U(393)2#1-0 
U(4_,4)=140 
U(595)2CT 
U(5e96)=-V¥*¥ST 
U(6e5)=V¥ST 
U(696)=CT 
RETURN 

END 


SUBROUTINE STIFOO 
DIMENSION U( 696) 


IF (RHO) lols2 
V=-1.9 

GO TO 3 

V=1.0 
CT=COSF(THETA) 
ST=SINF (THETA) 
DO 4&4 J=lo6 

DO 4& K=1+96 
Ut(Je«K)=0-0 
U(1lsl)=CT 

Ut ls4)=ST*V 
Ut(202)=CT 
U(295)=ST*V 
U(353)=1-0 

U( 491) 2=-ST*V 
U(454)=CT 
U(592)2-ST*V 
U(5e5)=CT 


(THE TAsRHOeU) 


(THETA sRHO oU ) 


010080 
010090 
010100 
010110 
010120 
010130 
010140 
010150 
010160 
010170 
010180 
010190 
010200 
010210 
010220 
010230 
010240 
010250 
010260 
010270 
010280 
010290 
010300 
010310 
010320 
010330 
010340 
010350 
010360 
010370 
010380 
010390 
010400 


‘010410 


010420 
010430 
010440 
010450 
010460 
010470 
010480 
010490 
010500 
010510 
010520 
010530 
010540 
010550 
010560 
010570 
010580 
010590 
010600 
010610 
010620 
010630 


10 


10 


U(6s6)=1-0 
RETURN 


END 


SUBROUTINE STAVEC(SVsNoBC) 
DIMENSION BC (6922) 


DO 10 J=1s6 
BC(JoN)=0-4 

IF CSV-ledlslos2 
BC(4e9N)=1.20 
BC(59N)=1-0 
BC(69N)=1.20 
RETURN 

IF (SV-22)39394 
BC(1lsN)=1.0 
BC(29sN)=1.-0 
BC(39N)=1.0 
RETURN 

IF (SV-32)59596 
BC(39N)=1-0 
BC(59N)=1.-0 
BC(69N)=1-0 
RETURN 

IF (SV-4e) 79798 
BC(1sN)=1.0 
BC(39N)=1.20 
BC(5eN)=1-0 
RETURN 

BC (29N)=1.0 
BC(39N)=1.0 

BC (69N)=1.20 
RETURN 

END 


SUBROUTINE STAVEO 
DIMENSION BC(6922) 


DO 10 L=1s6 
BC(L»N)=0.0 

IF (SV=1-0) 15192 
BC(29N)=1.-0 
BC(5eN)=1.0 
BC(60N)=1-0 
RETURN 

IF (SV=220)39394 
BC(1lseN)=1.0 
BC(39N)=1.20 
BC(49N)=1-0 
RETURN 

IF (SV~320)59596 
BC(20N)=1.0 
BC(4eN)=1.0 
BC(69N)=1.0 
RETURN 

IF (SV-40e0)79795 


10 | 


010640 
010650 


- 010660 


010670 
010680 
010690 
010.700 


_ 010710 


010.720 
010730 
010740 


- 010750 


010760 
010770 
010780 
010790 
010800 
010810 
010820 
010830 
010840 
010850 
010860 
010870 
010880 
010890 
010900 
010910 
010920 
010930 
010940 
010950 
010960 
010970 
010980 
010990 


-011000 


011010 
011020 
011030 
011040 
011050 
011060 
011070 
011080 
011090 
011100 
011110 
011120 
011130 
011140 
011150 
011160 
011170 
011180 
011190 


10 


20 


25 


}=1.-0 
BC(4eN)=1.-0 
}=1.0 


SUBROUTINE INVERT (AoNoDol oM) 


PROGRAM FOR FINDING THE INVERSE OF A NXN MATRIX 


DIMENSION A(25025) 9L(25)9M(25) 
SEARCH FOR LARGEST ELEMENT 
D=1.0 

DO8BO K=leN 

LIK) =K 

M(K)=K 

BIGA=A(KeXkK) 

DO20 [=KoeN 

DOZO J=KoeN 

IF (ABSF(BIGA)-ABSF(A( I 9J)))} 190920920 
BIGA=A(IeJ) 

LIK). =I 

M(K)2J 

CONTINUE 

INTERCHANGE ROWS 

J=L(K) 

IF(L(KI-K) 35935925 

D030 I=leN 

HOLD=-A(KolI) 

A(KeT)eAldel). 

A(JoI)=HOLD 

INTERCHANGE COLUMNS 

I=M(K) 

IF(M(K)—K) 45945,37 

DOGO J=loeN 

HOLD=-A(JoeK) 

A(JeoK) 2FAlJel1) 

A(JeI)*#HOLO 

DIVIDE COLUMN BY MINUS PIVOT 
DO55 [=leN 

IF( 1-K)50.55950 

A(T eK FACT eK ISI HAIK eK) ) 
CONTINUE 

REDUCE MATRIX 

DO65 T=leN 

D065 J=leN 

IFC I“K) 57565957 

IF(J-K) 60965960 
ACTsJIHACT oK) FA(Ke JI 4A TI oJ) 
CONTINUE 

DIVIDE ROW BY PIVOT 

DO7T5 J=leN 

IF (J-K) 70975970 
A(KesJ)=AI(K es JISAIK 9K) 
CONTINUE 

CONTINUED PRODUCT OF PIVOTS 


}02 


011200 
011-210 
011220 
011230 
011240 
011250 
011260 
011270 
011280 
011290 
011300 
011310 
011320 
011330 
011340 
011350 
011360 
011370 
011380 
011390 
011400 
011410 
011420 
011430 
011440 
011450 
011460 
011470 
011480 


-011490 


011500 
011510 
011520 
011530 
011540 
011550 
011560 
011570 
011586 
011590 
011600 
011610 
011620 
011630 
011640 
011650 
011660 
011670 
011680 
011690 
011700 
011710 
011720 
011730 
011740 
011750 


100 
103 
105 
110 
120 


125 


130 


150 


D=D#A(K 5K) 
REPLACE PIVOT BY RECIPROCAL: 
A(KeK)=leO/A(KoK) , 
CONTINUE 

FINAL ROW AND COLUMN INTERCHANGE 
K=N aii 
K=(K-1) 

LF(K) 1505150103 

T=L(K) 

IFC E-K) 12091208105 

DO110 J=1osN 

HOLD=A( JK) 

AlJoK)=-AlJel) 

A(J»o!)=HOLD 

J=M(K) 

IF (J-K) 10051009125 

DO130 I=1loeN 

HOLD=A(KolI) 

A(Keol)=~Al Jol) 

A(JoI)}=HOLD 

GO TO 100 

RETURN 

END 


SURBROUTINF HANGER (CLXsCLY9CT79U) 
DIMENSION U(606) 

DO 1 J=196 

DO 1 K=196 


“U(JsK)=040 


Ut(1s1)=140 
U(2s2)=140 
U(393)2=1-20 
U(454)=1-20 
U(595)=140 
U(696)=140 
U(493)=CTZ 
U(5—92)=-CLY 
U(6e1)=CLX 
RETURN 

END 


SUBROUTINE HANGEO (CTX sCLZ CTY 9U) 
DIMENSION U(696) 
DO 1 J=196 

DO 1 K=196 
U(J9K)=020 
Ut(1ls1)=1.0 
U(2s1)=CTX 
U(292)=1-0 
U(393)=1-0 
U(494)=140 
U(5e4)=CTY 
U(595)=1-20 
U(693)=-CLZ 


loys: 


011760 
011770 
011780 
011.790 


011800 


011810 
011820 
011830 


011840 


011850 
011860 
011870 
011880 
011890 
011900 
011910 
011920 
011930 
011940 
011950 
011960 
011970 
011980 
011990 
012000 
0%2010 
012020 
012030 
012040 
012050 
012060 
012070 
012080 


.012090 


012100 
012110 
012120 
012130 
012140 
012150 
012160 
012170 
012180 
012190 
012200 
012210 
012220 
012230 
012240 
012250 
012260 
012270 
012280 
012290 
012300 
012310 


Ta 


14 
15 


U(6»26)=1-0 
RETURN 
END 


SUBROUTINE CFIELO (AJT»RoZ sGeSORI sRHO»THETAs DoDI oFFACSA) 


DIMENSION A(606) 

PHI =THETA/Z 

RHOD=RHO 

IF (RHO-O02O0) llellsi2 

V=-12e0 

GO TO l 

V=1.0 

CP=COSF (PHI ) 

SP=SINF (PHI) 

W=1e0/( G #AJT) 
F1l=(PHI#CP+SP)/2-0 
F3=(SP-CP#PHI)/2.0 
F5=(2e0-2e¢0*CP-PHI*SP)/2-0 
F6=(2e0*#PHI+PHI #CP-320#SP)/2e0 
RHO=ABSF (RHO) 

DO 2 J=16 

DO 2 K=16 

A(JoK)=020 

IF (SORT 1135013014 

AREA=( 326141592654/4.0) *#(D**#2-DI * #2) #F FAC 


- SD=RHO# THETA/ (G#Z*#AREA ) 


GO TO 15 

S$0=0.0 

A(l»s1l)=CP 

Alls2)=(W#RHO#F] )—(RHO#F3#R) 
A(ls4)=V#SP - 
A(1s5)=V# (WHR) #* (RHO#PHI#SP)/2.0 
A(1s6) 2V#*(W+R) #(RHO#*#2 ) HF 3 
A(2e2)=CP 

Al295)=V*#SP 

A(206) 2V*RHO#(120-CP) 
A(301)2-A(296) 
A(302)=-12e0#A(1+6) 

A(303)=1.20 

A( 304) =RHO#SP 

A( 305) =((R#*( RHO##2 ) *PHI#SP)/2.0) -W#RHORFQZHFS 
A(3e6)=(R#*RHO##3#F 3) — (WHF 6*#RHOH*3) -—SD 
Al(Gel)=-V#SP 

A( 492) =-V#(WtR) *RHO#PHI *#SP/2 0 
A(4e4)=CP 
A(4e5)=(R#*RHO#F 1] )— (W#RHO*F 3) 
A(4e6)=A( 355) 

A(592)=-V#SP 

A(5e5)=CP 

A(596)3A( 394) 

A(656)=1.20 

RHO=RHOD 

RETURN 

END 


104. 


012320 
012330 
012340 
012350 
012360 
012370 
012380 
012390 
012400 
012410 
012420 
012430 
012440 
012450 
012460 
012470 
012480 
012490 
012500 
012510 
012520 
012530 
012540 
012550 
012560 
012570 
012580 
012590 


- 012600 


012610 
012620 


(012630 


012640 
012650 
012660 
012670 
012680 
012690 
012700 
012710 
012720 
012730 
012740 
012750 
012760 
012770 
012780 
012790 
012800 
012810 
012820 
012830 
012840 
012850 
012860 
012870 


_ 


10) 
ll 


SUBROUTINE POINO (AMUsAIX sATY sOMEGAsTL sZsSORI 53) 


DIMENSION B66) 
AM=AMU*¥TL/(2-1¢0) 

BDO 1 J=196 

DO 1 K=196 

B(JsK)=020 

BE(lsl)=1-0 
B(291)=-AM*(AIX*OMEGA) ##2 
B(292)=1-0 

B(393)=1-0 

B(4e4)=1-0 

IF (SORI-02)939810 
B(504)=-AM*(ATY*OMEGA) #*2 
GO TO il 

B(504)=0 

B(5s5)=1e0 

B(693) =AM*OMEGA##2 
B(6s6)=1-0 

RETURN 

END 


SUBROUTINE CFIELD (RoZsGsSDsRHOsTHETAsDsDI sE 9A) 
DIMENSION A(696) 


 PHI=THETA/Z 


RD=ABSF (RHO ) 
IF (RHO) 19192 


V=-1-0 
GO TO 3 
V=1-0 


CP=COSF (PHI ) 

SP=SINF (PHI) 

DO 4&4 J=196 

DO 4 K=196 

A(JsK)2=0.0 

A(ls1)=CP 

At(ls2)=-V¥*SP 
A(1ls3)=-V#¥RD*¥(120-CP) 

A( 194) =-V¥RD#* 2#R* (PHI-SP) 
AR=3e141592654* (D*¥*2-DI**2)/4.0 
W=1-0/(E*AR) 
F3=(SP-PHI*CP)*.5 
Fl=(PHI*CP+SP)*e5 
F5=(260-2e¢0*CP-PHI*SP) #25 
F6=(2e0*PHI+PHI *CP—-3 eO*SP ) *e5 
A( 125) =V*¥(RD*#WEPHI *#SP# e5-R*¥FS#RDE*S) 
A(ls6)=RD#W#F1+R¥F 6RRD**3 
A(2s1)=V*SP 

A(292)=CP 

A(293)2=RD*SP 
A(224)=R*RD#*¥2%(120-CP) 
A(295)=RD*#F3*(W+R*ERD# #2 ) 
A{226)=A(195) 

A(393)=1e9 

A(394)=RD#R*PHI 


|o5 


012880 


~ 012890 


012900 
012910 


- 012920 


012930 
012940 
012950 
012960 
012970 
012980 
012990 
013000 
013010 
013020 
013030 
013040 
013050 
013060 
013070 
013080 
013090 
013100 
013110 
013120 
013130 
013140 
013150 
013160 
013170 


013180 


013190 
013200 
013210 
013220 
013230 
013240 
013250 
013260 
013270 
013280 
013290 
013300 
013310 
013320 
013330 
013340 
013350 
013360 
013370 
013380 
013390 
013400 
013410 
013420 
013430 


NANAN 


END 


A(325)=A(294) 
A(306)=A( 194) 
A(494)=140 
A(4e5)=Al293) 
A(4e6)=A( 193) 
A(595)=CP 
A(5s6)=A( 192) 
A(695)=Al( 201) 
A(626)=CP 
RETURN 

END 


SUBROUTINE POINT (AMUsATY sOMEGAs TL oZoR1 5B) 


DIMENSION B(626) 
AMSAMU#TL/(Z2—-le) 
DO 1 J=1l6 

DO 1 K=196 
B(JeK)=040 
B(lsl)=1e0 
B(2e2)=120 
B(3e3)=1-0 
IF(RI)29293 
B(493)=-AM* (AIY*#OMEGA ) ##2 
B(494)=1-40 

B( 502) =AM*OMEGA##2 
B(5e5)=1-0 
B(691)=-B(592) 
B(696)=1¢0 

RETURN 


SUBROUTINE SFIELD (DI sTL»RoZ9GesSDeDsEoFFAC 2A) 


DIMENSION A(696) 
ARA=(34141592654/4 20) #(D¥*2-D] #*2) 
AREA=ARA#FFAC 

DL =TL/2 

DO 1 J=295 

A(1leJ)=0-0 

A(lel)}2=1e0 

A(1s6)=DL/(E*#ARA) 

A(291)=0e90 

A(2e2)=1e0 

A(203)2DL 

A(294)=DL**2*R/2.0 

IF (SD)29293 
A(295)=DL¥*3*#R/6e0-DL/ (GFAREA) 
GO TO 4 

A(295)2DL**3#R/6.0 

A(226)*0e0 

A(301)=020 

A(302)}=0.0 


106 


013440 
013450 
013460 
013470 
013480 
013490 
013500 
013510 
013520 
013530 
013540 
013550 
013560 
013570 
013580 
013590 
013600 
013610 
013620 
013630 
013640 
013650 
013660 
013670 
013680 
013690 
013700 
013710 
013720 
013730 
013740 
013750 
013760 
013770 
013780 
013790 


‘013800 


013810 
013820 
013830 
013840 
013850 
013860 
013870 
013880 
013890 
013900 
013910 
013920 
013930 
013940 
013950 
013960 
013970 
013980 
013990 





= > => 


A(293)21¢0 
A( 304) =2=DL#*R 
A(395)=A(294) 
A(396)=0e20 
DO 5 J=16 
A(49J)=020 
A(59J)=040 
A(69J)=000 
A( 494) 2140 
A(495)=DL 
A(595)7120 
A(696)21¢0 
RETURN 

END 


SUBROUTINE SFIELO (DI oTLeAJTsRoZ9GeSDRI sDeFFACSA) 
DIMENSION A(696) 
AREA=(30141592654/400) #(D¥#2-DI*H2) HFFAC 
DL=TL/Z 

DO 2 J=196 

DO 2 K=196 

A(JsK)2020 

A(lol1)=1-0 

A(202)2#1¢0 

-A(393)=120 

A(494)=1-20 

A(595)=120 

A(696)=1¢0 

IF (SDRI-Oo) 79798 

A( 396) =DL##3*R/600-DL/ (G*AREA) 
GO TO 1 

A( 396) =DL##*3#R/600 
A(1ls2)=DL/(G#AJIT) 

A(394)=DL 

A( 395) =DL¥#2#R/200 

A( 455) =DL#R 

A(496)2A(395) 

A(596)2DL 

RETURN 

END 


SUBROUTINE MATMUL (AsBoVU) 
DIMENSION A(696)9B(696) 9C (696) 9U(696) 
DO 1 J=196 

DO 1 K=196 

C(JoK)=000 

DO 1 L=196 

Ci JoK HCl IeK) +B (Jol) FA (LOK) 
DO 2 J=16 

DO 2 K=196 

U(JsK)=C(J9K) 

RETURN 

END 
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014000 
014010 
014020 
014030 


- 014040 


014050 
014060 
014070 
014080 
014090 
014100 
014110 
014120 
014130 
014140 
014150 
014160 
014170 
014180 
014190 
014200 
014210 
014220 
014230 
014240 
014250 
014260 


- 014270 


014280 
014290 


014300 


014310 
014320 
014330 
014340 
014350 
014360 
014370 
014380 
014390 
014400 
014410 
014420 
014430 
014440 
014450 
014460 
014470 
014480 
014490 
014500 
014510 
014520 
014530 
014540 
014550 


64 


67 


65 


64 


67 


65 


SUBROUTINE FINBRA(UsVVeVoeMAT 9A5Z) 


DIMENSION U(606) 9A(6 96) 9VV(693) 9V( 693) 


TYPE DOUBLE VVsV 

IF (MAT) 65964965 

DO 67 J=196 

DO 67 K=l¢3 

VideKkK)jy=05 

DO 67 L=1>+%6 

Vi JoK) =Vi Je K)FUCJoL) ¥VVILoK) 
RETURN 

N=Z-le 

DO 5 M=loeN 

DO 4 J=16 

DO & K=1l»93 

Vi JoK)=O65« 

DO & L=1l6 
V(JoK)FVitJeoKI4U(JoLIFVVIL 9K) 
DO 5 J=1s6 

DO 5 K=l»3 

VV(JeK)EV( JK) 

DO 7 J=1s6 

DO 7 K=1l93 

Vi JeoK)=O6 

DO 7 L=1s6 
VitJoKSVlUeKIFAl Jo LI *FVVIL OK) 
RETURN 

END 


SUBROUTINE FINMAT (UsUUsVoeMAT 5A0Z) 


DIMENSION U(626)s A(626)9V(6903) sUU(603) 


TYPE DOUBLE UUsV 

IF (MAT) 65964065 

DO 67 J=16 

DO 67 K=193 

V(JoK)=06 

DO 67 L=16 

Vi SJoK REV SoK)4U( Jol) FUU(L dK) 
RETURN 

N=Z—-le 

DO 5 M=loN 

DO 4 J=1s6 

DO 4 K2193 

Vi JoK) #0. 

DO 4 L=196 

V(JoK HVC JoK)4U( Jol) FUU(L 0K) 
DO 5 J=16 

DO 5 K=193 

UU( JS oK) SV Sok) 

DO 7 Jels6 

DO 7 K2=1l.e3 

Vi JeoK)=06 

DO 7 L=16 

VO SeKI FV SeKI+AlJoL) FUU(L 0K) 
RETURN 
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014560 
014570 
014580 
014590 
014600 
014610 
014620 


_. 014630 


014640 
014650 
014660 
014670 
014680 
014690 
014700 
014710 
014720 
014730 
014740 
014750 
014760 
014770 
014780 
014790 
014800 
014810 
014820 
014830 
014840 
014850 
014860 
014870 
014880 


014890 


014900 
014910 
014920 
014930 
014940 
014950 
014960 
014970 
014980 
014990 
015000 
015010 
015020 
015030 
015040 
015050 
015060 
015070 
015080 
015090 
015100 
015110 


END 


SUBROUTINE BRANCH (R39VVePHI su) 


DIMENSION R3(393)9VV(693) 9U(656) 9R1 (25925) s9LL (25) 9MM(25) 9R2( 393)» 


1R(393)2G61(393)9G2(3s3) 
TYPE DOUBLE vv 
DO 4 L=1+s3 
Rl(leLI=F=VVilsL). 
R1(2sL)=VV(2sL) 
R1(30L)=VV(39L) 
R2(leL)=VV(4oeL) 
R2(2esL)I=VV(5sL) 
R2(3eL)=VV(69L) 
SP=SINF (PHT ) 
CP=COSF (PHI) 
Gl(le1)=CP 
Gl(19s2)=-SP 
G1(193)=0. 
G1(2+%1)=SP 
Gl(202)=CP 
G1(293)=0. 
G1(391)206¢ 
G1(392)=06 
G1(393)=le 
G2(lel)=le 
G2(192)=06 
G2(103)=0-6 

-G2(291)=06 
G2(292)=CP 
G2(293)=-SP 
G2(391)=06 
G2(3%92)=SP 
G2(393)=CP 
CALL INVERT (R13sDsLL »MM) 
DO 5 J=ls3 
DO 5 K=1l9s3 
R3(J9K)=0.0 
DO & L=1] 3 
RZ(JoKI=ZFRZ(IUeK)FR1(JeoL)*GL(L Kk) 
DO 6 J=193 
DO 6 K=1l 93 
R1(JoK)=020 
DO 6 L=le3 
R1l1(JeKI=ERL(JeK)+G2(JoLI#FRZ(L 9K) 
DO 7 J=ls3 
DO 7 K=1+3 
R(JsK)=020 
DO 7 L=ls3 
R(JeKIHFR(JeoKItR1L (Jol) #RI(L 9K) 
DO 8 J=196 
DO 8 K=1+%6 
U(JsK)=020 
U(1ls1l)=1-0 
U(292)=1e0 
U(343)=1¢0 


015120 
015130 
015140 
015150 
015160 
015170 
015180 


° GLS5196 


015200 
015210 
015220 
015230 
015240 
015250 
015260 
015270 
015280 
015290 
015300 
015310 
015320 
015330 
015340 
015350 
015360 
015370 
015380 
015390 
015400 
015410 
015420 
015430 
015440 


(015450 


015460 
015470 
015480 
015490 
015500 
015510 
015520 
015530 
015540 
015550 
015560 
015570 
015580 
015590 
015600 
015610 
015620 
015630 
015640 
015650 
015660 
015670 


U(4eL)=R(1sL) 
U(5eL)=R(29L) 
Ul6esL)=R(3eL) 
RETURN 

END 


SUBROUTINE BRANCO (R3eVVePHI sU) 


DIMENSION R3(393) 0U(696) 9VV( 6593) sR1125925) 9R2(393)>5 


TYPE DOUBLE VV 


161(393) 9G2(393) 98( 393) sLL (25) »MM(25) 


DO 4 K=1l93 
R1(1lsK)=VV( 1K) 
R1(29K)=VV(3 9K) 
R1(39K)=VV(49K) 
R2( 19K )=VV(29K) 
R2(29K)=VV( 50K) 
R2(39K)=VV( 69K) 
SP=SINF (PHI) 
CP=COSF (PHI ) 
G1l(191)=CP 
~61(192)=0. 
G1(193)=-SP 
G1(251)=0.6 
Gll2s2)=le 
G1(293)=0.6 
G1(391)2SP 
G1(392)=06 
G1(393)=CP 
G2(191)=CP 
G2(192)=SP 
G2(153)=20. 
G2(291)=-SP 
G2(292)=CP 
G2(293)=06 
G2(301)=0.6 
G2(392)=0.4 
G2(393)=le 
CALL INVERT (Rls3sDsLt»MM) 
DO 5 J=193 
DO 5 K=1ls3 
R3(J0K)=0.20 
DO 5 L=l»3 
R3A(SesK)HR3A(JoK)+R1 (Jol) #*Gl(L ok) 
DO 6 J=l3 
DO .6 K=193 
R1(J9K)=0.0 
DO 6 L=1l»93 
R1I(JosKIRFR1I(UeK)+G2(JoL)IFRZ(LoK) 
DO 7 J=Als3 
DO 7 K=1l93 
S$(JsK)=020 
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015680 
015690 
015700 
015710 


- 015720 


015730 
015740 
015750 
015760 
015770 
015780 
015790 
015800 
015810 
015820 
015830 
015840 
015850 
015860 
015870 
015880 
015890 
015900 
015910 
015920 
015930 
015940 
015950 
015960 
015970 
015980 
015990 
016000 
016010 
016020 
016030 
016040 
016050 
016060 
016070 
016080 
016090 
016100 
016110 
016120 
016130 
016140 
016150 
016160 
016170 
016180 
016190 
016200 
016210 
016220 
016230 


WW 


fw 


DO 7 L=le3 


StU eK) SS(JoKIFRL(JoL)#RI(L 9X) 


DO 8 J=126 
DO 8 K=196 
Ut(JsK)=020 
Ut(lsel)=12e0 
U(291)=Stlel) 
U(292)21¢0 
U(293)=St1le2) 
U(294)=St193) 
U(393)=1-0 
U(494)=1.0 
U(591)=S(291) 
U(593)2=S(292) 
U(594)=S(293) 
U(5e5)=120 
U(6091)=St3e1) 
U(693)=S(392) 
U(694)=S(393) 
U(626)=1-0 
RETURN 

END 


SUBROUT INE STATEM(SVI sUU9D19029P19P2ePP»PMsSS) 
DIMENSION SVI (601) sUU(693) 
TYPE DOUBLE UUsP1sP29D1sD2 


M=0 

DO 1 J=16 

DO 1 K=1l»s3 
UU(J2K)=O06 

IF (PM) 40040941 
DO 42 J=1%6 


IF(SVI(J0l)) 42942943 


M=M+] 

UU(JeM)=Hle 

IF(M-1) 42944942 
UU(J02) 2-01 
UU(Js3)2-D2 
CONTINUE 

RETURN 

P1=P1+01 

P2=P2+D2 

IF (SS) 14915915 
IFtSS-le) 12911910 
IF(PP-1.e) 12911910 
DO 4&4 J=126 
IF(SVI(J91)) 49492 
M=M+] 

UUtJeM)=le. 
IF(M-2) 49395 
UUt(Je1l)=P1 

GO TO 4 

UUt(Je1l)=P2 
CONTINUE 


1 


016240 
016250 


~ 016260 


016270 
016280 


* 016290 


016300 
016310 
016320 
016330 
016340 
016350 
016360 
016370 
016380 
016390 
016400 
016410 
016420 
016430 
016440 
016450 
016460 
016470 
016480 
016490 
016500 
016510 


- 016520 


016530 
016540 


"016550 


016560 
016570 
016580 
016590 
016600 
016610 
016620 
016630 
016640 
016650 
016660 
016670 
016680 
016690 
016700 
016710 
016720 
016730 
016740 
016750 
016760 
016770 
016780 


POTLIRN 
fee 2) J=)]95 
BetoVitJs1)) 2ls2ls22 


‘s= M+) 

UU(Js4%)=le 

TF (M=-1) 25926525 
UU(Jo3)=P1 
PE(M—7) 219243921 
HEC Js 2) =O 
CONTINUF 

RETURN 


DO el J=135 

IF (SVI(J91)) 31931932 
M=M+] 

YUCIJIsM)=1. 

IF (M-1) 35536935 
UU(JI22)=P2 
[F(M=-7) 2153] 233 
UU(Js2)=P) 
CONTINUE 

RETURN 

ENT 


SUBROUTINE STATFRI SVR oe VV oN os YYV ok sp BARC oD] 92 PP 9 P49 $5501 502903 5F INK) 


DIMENSION VV(693) »¥Y(3939225) .V(393) sSVB(5222) 
TYPF DOUBLE vv 
31=Ge 
B2=0. 
B3=06« 
Pi= ) e 
De =O'5 
IF (PM) 32394%922 
2 TF( PRO) 21470921 

meet) } 73570922 
Hy IO” sal > 
eel) K=1,52 
VO JesK)HVYVY(UoK oN) 
Pet ss) Lllsl2s1l2 
IF(PP-1.) 32931930 
IF (SS—-le) 33931559 
BIHVOlsl)+VC192)*DO14V01 52) #0248) 
B2=V(291)4V(292)%*D14V (203) *D24+R2 
RI=V(AgTP4V( 397) EDLFEV( 252) FN74R 
Beet 42 
RIHRL HAVE 1,1) #AN24YV (1,7) #1 +N (1,7) 
R2=FR24V (2 6)-) #024 (7,72) #N14V(2 2) 
PI=R34V(451)*#D24+V( 7,2) 8#D14+V( 2,2) 
Gor {TO 22 
BLl=Sl4V( 1-18 014+V( i922) ¥O24+V(153) 
F2=024V(291)*D14V( 292) #024+V (252) 
GP o=-34V(0501)*514V( 352) *024+V( 452) 
el =) baer) . 
Ween lrOes+e2 
VAS QAlLFDN354N 2 


se 


is 


016799 
016800 
016810 
016820 
016830 
016840 
016850 
016860 
016870 
16880 
916899 
016900 
916910 
016920 
916930 
016940 
016950 
016960 
016970 
016980 
916990 
A1LTANA 
917019 
917929 
917039 
017040 
017050 
017060 
017070 
017069 
917090 
917199 
O17119 
917129 
017139 
017149 
017150 
017169 
Oe? 1 70 
017180 
017199 
017200 
One? 210 
G1 7220 
917229 
917249 
Ol 250 
017269 
017270 
017280 
A790 
017309 
O21 7310 
917329 
£17330 


20 


51 
41 
43 


44 


42 


20 
10 
12 


11 


M=1 

DO 1 J=1s6 

DO 1 K=1l93 
VV(J9K)=0.0 

DO 4 J=1s6 
IF(SVB(JeN)) 49492 
IF (M—-1)79397 

IF (PM) 345935935 
VV(J01)=01 

GO TO 36 
VV(Jeoll=le 

M=M+1 

GO TO 4 

IF (M=2) 49899 
VV(J01)=Q2 
VV(JoM)=ele 

M=M+1 

GO TO 4 
VV(J01)=Q3 
VV(JSJoM)=le 
CONTINUE 

RETURN 

IF(FINK) 54951954 
IF (BRC) 51951952 
DO 53 J=1le3 

DO 53 K=le3 


r VISeoKIVFYV (Je K oN) 


DI=(VELe2ISViLlolIFVE2sZI/Vi2e1l)4Vi3e2d/Vi 301))/36¢ 
D2=(V(193)/V0191)4V0293)/V(201)4V(393)/V(391))/36 
M=0 

DO 41 J=196 

DO 41 K=193 

VV (JSeK)=06¢ 

DO 42 J=19s6 

IF(SVB(JeN)) 42942943 

M=M+1 

VV(JoM)=le 

IF (M—-1) 42944942 

VV(J02)=-Dl 

VV(J93)==-D2 

CONTINUE 

RETURN 

END 


SUBROUTINE BRCOR(R3sUUsXXeJKeVVePM) 
DIMENSION R3(393) 9UU(6 93) 9XX(393) sUUU( 393) sVV( 693) 
TYPE DOUBLE VV.».UU 

IF (PM) 20921920 

IF( JK) 10910911 

DO 12 J=1ls3 

DO 12 K=1l93 

UUU (JsK)=UU(IsK) 

GO TO 3 

DO 13 K=l93 

UUU(1sK)=UU(19K) 


hls 


017340 
017350 


017360 


017370 
017380 


* 017390 


017400 
017410 
017420 
017430 
017440 
017450 
017460 
017470 
017480 


017490 — 


017500. 
017510 
017520 
017530 
017540 
017550 
017560 
017570 
017580 
017590 
017600 
017610 
017620 
017630 


017640 


017650 
017660 
017670 
017680 
017690 
017700 
017710 
017720 
0177390 
017740 
017750 
017760 
017770 
017780 
017790 
017800 
017810 
017820 
017830 
017840 
017850 
017860 
017870 
017880 
017890 


UUU(29K)=UU(3» 
UUU( 39K) =UU( 4» 


DO 1 J=1»s3 

DO 1 K=193 
XX(J9K) #04 

DO 1 L#1l»s3 
XXCSoK) XX UeoK) tR3( Job) #UUUIL»K) 
RETURN : 
IF (SUK) 25525922 
DO 23 J=1e3 

DO 23 K=19s3 
XX(JoK EVV SeK) 
RETURN 

DO 24 K=1l93 
XX(1lsK)=VVC( 1K) 
XX(29K) 2VV(39K) 
XX(30K)=VV(49K) 
RETURN 

END 


K) 
K) 


SUBROUTINE DELMA(UUs SVE sV9X29Y2sKKK 9X1 FF »PP sHH»sRRoSSeV1l sREMsFINK» 
1PMs01sD02) 


DIMENSION UU( 693) 9V(693) »SVE(651) 90(393) 
TYPE DOUBLE UUsV9D9X1sX2sY1sY290V901 902 
L=0 

X1=0-6 


 X2206 


111 


Y2=0-6 

Y1=06 

01206 

D2=0. . 

DO 111 J=13 

DO 111 K=1ls3 
D(JoK)=06 

DO 91 N#1 6 
IF(SVE(Ns1)) 9192991 
L=l+1 

DO 92 K=1s3 

Vib »sK)=UU(NsK) 
CONTINUE 

IF(PM) 665964561 
IF(FINK) 63962963 
O1=(V(192)/V0191)4V0292)/V(291)4+V(392)/V(391))/36 
O02=(V(193)/V(191)4V(293)/V(291)4V(393)/V(391) 30 
FINK=l. 

RETURN 

FINK=0.¢ 

RETURN 

RETURN 

IF (SS) 311293122312 

IF (PP-1.) 40024015402 
IF(SS—le) 31393149315 
IF(RR~le) 30029403300 
IF(RR-1le) 34053415340 
IF(V(393)) 93994593 


114 


017900 
017910 


017920 
017930 
017940 
017950 


' 017960 


017970 
017980 
017990 
018000 
018010 
018020 
018030 
018040 
018050 
018060 
018070 
018080 
018090 
018100 
018110 
018120 
018130 
018140 
018150 
018160 
018170 
018180 
018190 
018200 


-018210 


018220 
018230 
018240 
018250 
018260 
018270 
018280 
018290 
018300 
018310 
018320 
018330 
018340 
018350 
018360 
018370 
018380 
018390 
018400 
018410 
018420 
018430 
018440 
018450 


340 
93 


95 


94 


403 
341 
96 


97 


98 


102 


101 


100 
354 


355 


IF (V(392)) 93994593 

be.95 J=293 

DO 95 K=1.3 

D(JeK)2VJ5eK) 
DV=V(1902)#D(293)-V(193)*D(292) 
IF(DV) 98594598 

RR=1¢e0 

HH=le 

RETURN 

IF(V(293)) 969200996 

IF (V(202))9639200596 

DO 97 K=1e3 

D(2eK)=V(39K) 

D( 3K) =V(25K) 
DV=V(102)*D(293)-V(193)*D( 292) 

IF (DV) 985200998 
X1=X1—(V0191)*D(253)-Vl193)*#D(251))/DV 
X2=X2~-(V0192)*D(291)-Vl1l91)*D(292))/DV 
IF (X1) 10051029100 

IF(X2) 10091019100 

FF=1-0 

RETURN 

IF(SS) 35453559355 
Y2=Y2-D(391)/D( 393 )—-X1*D(392)/N(393) 
RETURN 
Y1=Y1~-D(391)/D(392)-X2*D(393)/D(392) 


. RETURN 


99 
20 


aT 
23 


314 
401 
330 
301 

81 


83 


82 


coa 
404 
84 


85 


86 


IF(SS) 20921921 
PP=l. 
GO TO 23 
SS=le 
RR=06 
HH=l. 
RETURN 
IF(RR~-le) 33093315330 
IF (RR~-120) 30194045301 
IF(V(391)) 815982581 
IF(V(392)) 812982581 
DO 83 J=293 
DO 83 K=1l93 
D(JoK)=VE eK) 
DV=V(101)#D(202)-Vlle2)*D( 291) 
IF(DV) 86282586 
RR=1¢e0 
HHele 
RETURN 
IF(V(291)) 84970984 
IF(V(202)) 84270984 
DO 85 K=1l93 
D(29K)2V(390K) 
D(3eK)=V(29K) 
DV=V(101)*D(292)-Vl192)*D(291) 
IF(DV) 862970586 
X2=X1—(V0193)4#D(292)-Vl1922)*D( 293) )/DV 
X2=X2—(V(191)*#D(293)-Vl1l93)#D( 291) )/DV 
IF (X1)2039201 9203 


ry ger 


HIS 


018460 
018470 
018480 
018490 
018500 
018510 
018520 
018530 
018540 
018550 
018560 
018570 
018580 
018590 
018600 
018610 - 
018620. 
018630 
018640 
018650 
018660 
018670 
018680 
018690 
018700 
018710 
018720 
018730 
018740 
018750 
018760 


(018770 


018780 
018790 
018800 
018810 
018820 
018830 
018840 
018850 
018840 
018870 
018880 
018890 
018900 
018910 
018920 
018930 
018940 
018950 
018960 
018970 
018980 
018990 
019000 
019010 


313 


405° 


211 
212 


213 
351 


350 


IF(X2) 20392029203 

FF=1.0 

RETURN 

IF (SS) 3523539353 
¥1=Y1-D(393)/D(391)-X2*D(392)/D(391) 
RETURN 
Y¥22=Y¥2-D(393)/D(392)-X1*D(391)/D00392) 
RETURN : 

IF(SS) 25924924 

SS=2e 

GO TO 26 

PP=2. 

RR=O6« 

HH=le 

RETURN 

IF(RR-le) 32093219320 

IF(RR=-le) 30294059392 
IF(V(303)) 7lse72s71 
IF(V(391))71972971 

DO 73 J=293 

DO 73 K=1l¢93 

DC JoK)=VC eK) 
DV=V(1e1)#D0(203)-V(123)*#D(291) 
IF(OV) 76972976 

RR=1.0 

HH=1le 

RETURN 

IF(V(201)) 74999974 

IF(V(293)) 74999974 

KKK=] 

RR=0-6 

RETURN 

DO 75 K=l»93 

D( 29K) 2V(30K) 

D(3oK)#Vl20K) 
DV2zV(191)*0(293)-V(1ls3)*D(201) 
IF (DV) 76999976 
KX2=X2-(V0192)*D(203)-V(193)*D( 292) )/DOV 
X1=X1-(V101)*0(292)-D(201)*V(102) )/DV 
IF(X1) 21392115213 

IF(X2) 21392129213 

FF=1.0 

RETURN 

IF(SS) 350+351.351 
Y1=Y1-(0(392)4X2*D(301))/D0(393) 
RETURN 
Y2"Y2-(0(352)4+X1#*D0(393))/D0(391) 
RETURN 

END 
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019020 
019030 
019040 
019050 
019060 
019070 
019a80 
019090 
019100 
019110 
019120 
019130 
019140 
019150 
019160 
019170 
019180 
019190 
019200 
019210 
019220 
019230 
019240 
019250 
019260 
019270 
019280 
019290 
019300 
019310 
019320 
019330 
019340 
019350 
019360 
019370 
019380 
019390 
019400 
019410 
019420 
019430 
019440 
019450 
019460 
019470 
019480 
019490 
019500 


APPENDIX E 
ACCURACY OF METHODS & SAMPLE PROBLEMS 
E-1 General remarks. 

Appendix E contains the results of vibration analyses of several 
piping systems. The first group was performed to establish the accuracy 
of the methods and compare them with the VIPIPE output. The second group 
was performed to establish assurance of solution by comparing the natural 


frequencies got from various methods. 


E-2 Method accuracy analyses. 
E-2-1 Straight sections. 
End conditions. 
System 1 - fixed-fixed. 
System 2 - fixed-free. 
System 3 - fixed-propped. 
System parameters (for system 1,2,3). 
length - 200 inches. 
diameter - 2.0 inches. 
wall thickness - 0.125 inches. 
Intensive properties (for system 1,2,3) 
density - 490 lbs/cu-ft 
shear modulus - 12x10° psi 
elastic modulus - 30x10° psi 
Mathematical model: Distributed mass with the effect of shear 
deflection and rotational inertia neglected to permit comparison with 


classical solution. 
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Comparison of results. 

As explained in Appendix B-9-4, when the iteration reaches the 
solution acceptability criterion, the corresponding natural frequency 
is determined using linear interpolation. Theoretically exact compari- 
son values are obtained by use of the same iteration method with the 
solution acceptability criterion as 0.0007 rad/sec. The values obtained 
using VIPIPE have a solution acceptability criterion 0.0185 rad/sec. 
Comparison is made by percent difference. 

The M-D method fails the systems 1 through 3, because the purging 
factors turn out to be all zeros. So for the purpose of testing the 
M-D method, modified System 1 is provided below. 

The letters D or S in parentheses denote double precision and single 


preaision respectively. 


Modified System 1: 

Add a very small curved section to the System 1. Second section is 
a mathematical model which has the same property as System 1 except the 
system parameters. 


Second section parameters: 
ve 


; / included angle of arc - 45 degrees 
f radius of curvature - 0.01 inches 
(This value of radius of curvature is 
not possible in an actual system, since the radius of curvature less than 


the radius of the pipe. In effect subroutine STIFCO or STIFOO will pro- 


vide the transfer matrix of a massless rigid corner). 
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Results and comparison.* 


a. System 1 and modified System 1.** (in-plane). 

































Mode |: | —~—=—=—~—s*Freq.(rad/sec) and difference in% sd 
No. Comparison value Delta m. (D) P-M. -A(D P-Mm. -F(D) | 
pol 7509994540 75.09994544 75.09994550 75.09994571 | 
| | 0 z 0 0 4 
ine 2 207 .01589036 |  207.01589094 207 .01589228 | 
| 5.0x107/% | 0 z 4.4x1077 2 
| 3 405.83391916 = | 405. sill | 405 83403022 405.83394702 | 
2 | enor, 2.4x10° 7% 7.0x107'% 
4 | 670.86408383 | 6 670.86310856 67086341454 67086345124 
| | | 2.8x10°% 1.0x1077%, 9.1x107°% | 


ag en ee -————— 
Ts | 1002.15520225 | -:1002.19978814 =| —- 1002. 18735784 1002.19300449 
:) | 4.7x107 74, 3.2x10 a 3.8x10°-%, | 
| 6 1399 .70436269 | 1399.59981662 | 1399, 64448571 


{ : a | - 

| : 7.5x107~% 4.3107 7% 
195959556134 

5.1 % 




















et: ee. 4 oe S te eee ig A =~ 


7 : 1863.51172599 ) 
| | 








ry 





ce ee eee ee EN ne 


 M-D m. (S)*** 
75. 09994536 | 










- -75.09994865 


| 75.09994539 
{ °, ‘ os °, 
me | oO *% | 4.0x10 % 




































<i Se ee 
2 | 207.01589228 | -207.01589071 +—=—S—=«207.01588836 |  —207.01588891 
| 5.0x10°°% 0 a 1.5x10° % | 0 z, 
(3 | 405. 83398145 ~-405.83403455 | 405.83395151 | 405.83401120 | 
i | 1.5x10° a | 2.1x107% | 9.5x10°° % | 2.0x10°? % 
4 | 670.86345124 | 670.86256288 |  670.86374475 |  670.86118054_ 
7.5x10°°% 2.2x10°* % | 5.1x10? % > 4.3x107" % 
“5 | | 1902. 23499946 | 1002. 19909306 1002. 21779856 
‘4 oe | .0x107> 4 — 4.4x107" % 6.3x10"> % 
| 6 | 1399. 17084415 1399, 94044241 
_ | asx"? 2 | eto 


*Sources of comparison frequencies for system 1 through 3 are given in E-2-2. 
(Additional footnotes on next page). 
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**In all the tables which follow, simple and evident abbregjations are used. 
For example, P-M m.-A(D) means the P-M method, sub-procedure A in double 
precision arithmetic. 


*kkModified System l. 
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b. System 2.(in-plane) 


Mode ws 
NO Comparison values 
1 11.80 213586 





73.96272297 
- 207.09776595 
405.82896584 
|” 670.86435913 
1002.1551877 
1399.70436329 


| 8 | 1863.51172626 
{ 


c. System 3. (in-plane) 














_Freq. (rad/sec) and difference in 2 


ee 


Delta m. (D) 


a ee 


11. 80213389 
5.5x107? % 
73.96272167 
1.7x10~© ¢ 
207 .09776685 
4.1x10-9 % 
405.82894132 
6.1x10~° % 


670 .86200817 
1002. 26833752 


1399 .08849534 





-2,.0x107- % 





186 2.559837 28 


P-M m.—A (D) 


A oe 


11.80214614 
Je 5x10"° — 
713.9627 2478 
_l.Tx10-6 # 
207 .09776388 
ee ne SALA 
405.82893748 
7.0x10"° % 
670 .86326697 
1.7105? Ae 





100 2.18809484 


4.4x1072 % 


1399.61231065 


6.6x107? & 







5.1x107° % 


ee __ Fens. (rai/ons) and difference in % 
NO Comparison values Delta m. (D) P-M m.-A (D) 





P| 51.75397 285 


167.71502090 
349..92609071 
4 598. 39432089 


913.12074582 
Z& 1294.10536551 
1741.34817946 








we ee 


51.75397263 
0 ia 

167.71602009 
Oo. 

349 .92608231 
1.1x107/ % 

598. 39381063 


8.5x107? % | 


913.13150640 
1.2x107? % 
1293.77838653 


2.75 10B soa at 


51.75397569 

0 %o 

167.71602200 
O % 

349 .92610259 
_3.4x10~° 

598. 39446341 
2.5x1079 % 


~ 913.13180993 


7 1. . 2x10"? % 
1 294, 11223844 
__5.8x107* % 
1748. 59011823 


4.2x107! & 


— 





—s a 


PM m.-F (D) | 
11.80213746 
1.4x107? 4 
73.96272367 
0 qe 
207.09776562 
0 & 
405.82900191 
g.1x10~° % 
670.86325735 
1.7x10-9 % 
1002.19715557 
4,2x107-2 &- 
1399.14599562 
4.2x1077 & 







a 
’ 


P-M m.-F (D) 
51.75397564 
0 % 
167.71602189 
Oo 
34992610209 
3.4x10-© % 
598 . 39407344 
4.2x107? % 
913.13109335 
1.2x1077 % 
1294.11761177 
1.9x107* % 


E-2-2 Sources of comparison frequencies for system 1 through 3. 

The governing fourth order equation of a homogeneous single component 
straight system , 

so AW 
was solved to give 
y = a cos dx + b coshXx + ¢ sin dA*' t+ 4 sinhXx 

Upon applying boundary conditions for System 1 (y = y' = 0: atx =0, 
x = L) and simplifying, it was found that the eigenvalues of the system 
must satisfy 


cos 0: cosha = 1 (E-2-2-1) 


where -- +f wt a 
EL (E-2-2-2) 


Upon applying the boundary conditions for System 2(y = y' = 0 at x = 0, 





y" = y" = 0 at x = L), it was found that the eigenvalues of the-system 
must satisfy 

cos + coshX = -1 (E-2-2-3) 
Upon applying the boundary conditions of System 3(y = y' = 0 at x = 0, 
y = y" = 0 at x = L), it was found that the eigenvalues must satisfy 

sind - coshX- cos &% :- cosh & = 0 (E-2-2-4) 
The roots of equation E-2-2-1, E-2-2-3 and E-2-2-4, from which the natural 
frequencies of the corresponding systems may easily be found, were ob- 


tained through iteration using a digital ‘computer. 
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E-2-3 Curved section (System 4) 


End conditions: 


Mathematical model; 


Results. 


a. In-plane vibration* 






369.11199342 











571.22235666 


System parameters: 


50 inches 
270 degrees 

2.0 inches 
.125 inches 


radius of curvature 
included angle of arc 
diameter 

wall thickness 


Intensive properties: 


density 556. lbs/cu-ft 
shear modulus 12x106 psi 
elastic modulus 30x10® psi 


fixed-fixed. 
lumped mass with shear deflection and rota- 


tional inertia neglected. 











_ Mode . Frequency (rad/sec ; l 
| No. M-D m. (D) P-M m.-A (D)| P-Mm.-C (D) © 
1 «| 66.69934902 66.69934902 | 66.69934980 66.69934948 
| 2  1171.22352538 171. 22352538 | 171.22352610 | 171.22352592 
4 3  |335.35007083 | 335.35007083 335.35007083 | 335.35007174 
mi 535.64529414 | 535.6429665 | 535.64529637 
: 5 773.44358487 | 773.44358842 | 773.44358812 
b. Out-of-plane vibration* " 
{Mode | a _Frequency (rad/sec) | 
| No. | Delta m.(D) M-D m. (D) | P-M m.-A (D) | 
3536309735 35.36309735 | 35.36309778 
: | 94.42588921 | 94.42588921 | 94.42588653 - 
: 207.89287257 | 207.89278647 | 207.89287345 _ 
nnn <— 
Pd 369.11278647 | 369.11278689 
; - — aan 
| 


preven i 
wr 


*See footnote next page. 


571.22235577 
ee 
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*The fundamental frequency value of in-plane and out-of-plane case are, 
close to the values obtained by employing equation W= H(o):( EI/uR )” 
after obtaining F( « ) from the curves in reference 


[5]. Natural 
frequency values from the three methods are the same up to six digits. 


One can get more higher mode frequencies using the M-D and P-M method. 
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E-3 Other Sample Problems 
A. System 5 


System parameters 


a tle ee 







a component - 
| Parameter 





“dianeter ( (in) |2-0 
A 
‘wall thickness .).75 0.75 0.75 N 0.75 0.75 0.75 
| (in) ! G 
‘length (in) ‘130. 12g) 620 E 60.21 1.50 6.28 
R 
radius of curve | 4.0 30. 4.0 


_ Cin) : 
| inel. angle of are| 20. 115. 90. 
i ee ee a 


Intensive properties of all components except hangers: 





density: 556 lbs/cu-ft.* 
shear modulus: 12x10° psi 
elastic modulus: 30x10° psi 
Hanger spring rates: 
CLX , CLY ,CLZ 10001bs/in. 
CTX , CTY ,CTZ 1000 in-1b/rad. 
End conditions: 
both ends of system fixed. 
Mathematical model: Distributed mass with the effects of shear deflec- 


tion and rotational inertia considered. 


6 f 








Nur 


} + 


~— 
5) er atl 


*This fictitious valye was chosen to permit comparing with results of 
Fink [4] who used the same value, 
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Results. (in-plane vibration) 



























7 Frequency(rad/sec) 
(a ee a a a i aaa: a lh 
Delta m. (D) N-D Me (D) P-M m.-A (D) P-M m,-C (2) 


116.89479180 


ne ee ee a 


2 304 .88354398 


i 3 566.43670976 


} mi 


mG 1666.58935151 


826 .882819 29 


— ———— - ee <e 


1126 .83797508 





eee ae el 


1126 .83799669 


1666. 58812881 


ee a= = eee oe en te er ee 


2099 .5352509 


ee 


a er ee ee oe 


1666. 58878192 


2099. pallies : 
: 








7 | | 2099. 54173011 
eee. 2 eae 





‘Results in single precision arithmetic and comparison with double precisbon results. 


Frequency(rad/sec) and diff, in % 


Mode wae nS 
NO Delta Mm. (s) | M-D nm. (8) P-M m.-A (S) 


P-M m.-C (S) 
























1 116.89479181 ; 
' Oo fF 
488655466081 (°° ~S—=<“<‘<i<S SCS 
1.0x107° % | | i 7 
5 56643620786 ‘on | 
8.9x107> % 2 a 
4 826.88281494 |  826.88281417 , 
1.0x107© % 6.0x10-7 1.9x10-7 4 
1126 .83768564 1126 .83786333 112683745223 
= = 4.4x10=° % 1.2x107? % 4.6x107> % 
1666. 58927426 | | 
| 7x10~° fo 
i Meco. Up | ia we 
5. ox10~ % > 


ee ee ee ee. 


—— merece Se ae ec = J _ 
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oe ee 


B. System 6 (single branch system) 


System parameters; 





components 





parameter 






dia. (inch) 







4- wall thickness 


| (inch) 
4—* 3° 2 a2 ‘length (inch) 


fi 
Branch 
incl. angle of 


arc. (degree) 
Intensive properties: 
density 556 lbs/cu-ft 


shear modulus 12x10° psi 


elastic 30x10° psi 
modulus 


End conditions: All ends fixed. 
Mathematical model: Distributed mass with the effects of shear deflec- 


tion and rotational inertia neglected. 


Results | 
Mode | _ Frequency (rad/sec) 
No. ; Delta m.(D) | M-Dm. (D) [| P-Mm.-A (D) | P-M m.-C (D) 

1 19425470558 | 
2 278.61078241 | 
3 628.59556355 | 
4 747.68138190 747 .68087208 | 747.68137310 74768137759 
5 1300.29478437 | 1300.31646085 1300. 31641281 
6 1381.79339427 | 1381.86234251 | 1381.86230734 
7 1991.80356348 1991.38946110  1991.38946101 
8 


2257.92926800 , 2257.58161443  2257.58016324 


a ee ee pee | 


jar 


| 
| 
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